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ABSTRACT 

A  technique  for  the  use  of  semi-infinite  planes  to 
approximate  the  log  magnitude  versus  log  frequency  charac- 
teristics of  two-dimensional  filters  is  presented  (Extended 
Bode  approach) .   This  technique  is  applied  to  quarter  plane 
filters  and  works  well  for  separable  transfer  functions. 
Other  non-separable  canonic  (basic)  transfer  functions  are 
also  studied  in  terms  of  their  planar  approximation.   The 
general  technique  is  shown  to  be  useful  for  the  insight  it 
provides  as  well  as  being  a  simple  approach  to  design.   The 
double  bilinear  z-transform  is  studied  for  use  with  the 
two-dimensional  analog  transfer  functions  to  convert  them 
to  the  digital  domain. 
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I.   INTRODUCTION 

A.   DIGITAL  SIGNAL  PROCESSING 

Digital  signal  processing  has  its  roots  in  16th  century 
mathematics,  especially  in  the  fields  of  astronomy  and  the 
compilation  of  mathematical  tables.   Today  it  has  become 
a  powerful  tool  in  a  multitude  of  diverse  fields  of  science 
and  technology.   The  applications  of  digital  signal  pro- 
cessing varies  from  low-frequency  spectrum  seismology 
through  spectral  analysis  of  speech  and  sonar  into  the 
video  spectrum  of  radar  systems  [1] . 

In  their  book,  Theory  and  Applications  of  Digital 
Signal  Processing  [2]  ,  Rabiner  and  Gold  have  combined  digi- 
tal signal  processing  theory,  with  a  variety  of  applications 
ranging  from  sonar,  radar,  communication,  music,  seismic 
and  medical  signal  processing,  with  digital  component 
technology.   This  technology  is  the  main  driving  force  for 
progress  in  this  field  as  well  as  the  general  area  of  com- 
puter design.   It  was  also  observed  by  them,  that,  although 
the  formulation  of  engineering  problems  is  often  as  vague 
as  those  of  the  "softer"  sciences  (such  as  anthropology, 
psychology,  etc.),  the  execution  of  these  problems  appears 
to  depend  on  greater  and  greater  accuracy  and  reproducibility 
The  capability  of  digital  systems  to  achieve  a  guaranteed 
accuracy  and  essentially  perfect  reproducibility  is  very 
appealing  to  engineers. 


A  digital  filter  is  defined  in  [3]  to  be  a  computa- 
tional process  in  which  a  sequence  of  numbers  acting  as 
input,  is  transformed  into  a  second  sequence  of  numbers, 
or  output  digital  signal;  where  the  term  "digital"  implies 
that  both  time  (the  independent  variable)  and  amplitude  are 
quantized. 

The  field  of  one-dimensional  digital  filtering,  which 
is  outlined  in  figure  1.1,  encompasses  recursive,  non- 
recursive  and  Fast-Fourier  Transform  processing.   The  terms 
recursive  and  nonrecursive,  instead  of  IIR  (Infinite  Impulse 
Response)  and  FIR  (Finite  Impulse  Response)  are  preferred 
in  this  thesis.   It  is  shown  that  recursive  processing  is 
much  more  efficient  than  nonrecursive  processing.   Stockham's 
method  [4]  to  perform  fast  convolution,  which  later  became 
known  as  the  FFT  method,  improved  the  efficiency  of  non- 
recursive  techniques,  so  that  comparisons  in  one-dimension 
are  no  longer  strongly  biased  toward  recursive  methods  [2] . 

One-dimensional  recursive  digital  filter  design  depends 
strongly  on  the  effective  and  well  developed  continuous 
filter  design  theory.   Moreover,  the  stability  analysis  of 
higher  order  one-dimensional  recursive  realizations  can  be 
solved  by  investigation  of  the  root  distribution  of  its 
factored  form  representation.   It  is  important  to  realize 
that  continuous  domain  design  techniques  and  factorization 
property  (the  fundamental  theorem  of  algebra)  exists  only 
in  one-dimension. 
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Areas  of  application  of  one-dimensional  digital  signal 
processing  are  listed  in  figure  1.1.   There  are  important 
military  and  specifically  naval  applications  in  addition 
to  those  listed,  such  as  missile  and  torpedo  guidance, 
launcher  control  system,  fire  control,  combat  information 
collection,  dissemination  and  display. 

Digital  filtering  in  several  dimensions,  which  is  over- 
viewed  in  figure  1.2,  has  gained  considerable  importance, 
especially  for  the  two-dimensional  case.   Figure  1.3  presents 
an  example  of  two-dimensional  low-pass  and  high-pass  filtering 
Until  196  6,  two-dimensional  digital  filtering  was  implemented 
by  nonrecursive  or  convolutional  techniques.   In  this  method 
the  output  is  the  weighted  sum  of  unit  sample  responses  of 
all  past  input  values.   The  serious  disadvantage  of  the 
convolutional  method  is  the  requirements  of  a  very  large 
number  of  arithmetic  operations. 

The  development  of  the  Fast  Fourier  Transform  (FFT)  in 
1966,  reduced  the  number  of  arithmetic  operations  consider- 
ably and  is  used  extensively  today.   Filtering  via  FFT  is 
accomplished  by  computing  the  transform  of  the  input  func- 
tion, multiplying  by  the  frequency  response  of  the  filter, 
and  inverse  transforming  the  result.   The  recursive 
algorithm,  which  has  in  general  an  infinite  impulse  response 
[3] ,  constitutes  another  technique  for  the  realization  of 
two-dimensional  digital  filters. 

Hall  [5]  compares  the  amount  of  computation  required  for 
the  three  filtering  techniques.   He  shows  that,  in  general, 
the  FFT  and  recursive  algorithms  are  preferable  to  the 

11 


02 

T 

*—* 

co 

o 

S 

1-1 

s 

CO 

En 

3 

en 

W 

>H 

2 

co 

H 

Q 

J 

1 

< 

Z 

Eh 

H 

s 

O 

H 

M 

Q 

•r1 

K 

fa 

O 

O 

fa 

CO 

>H 

2 

05 
O 

£ 

fa 

Eh 

K 

1 

Eh 

CS3 

i 

1 

fa 

OS 
D 
Eh 
U 
D 
« 
Eh 
CO 


CO 
H 
CO 

>• 
< 

z 
< 

fa 

Eh 

< 
Eh 
CO 


05 

a 

ij 

H 

&b 

§£ 

H  CO 

a 

z; 

o 

z 

eg 
fa 

Eh 
CO 

>H 

CO 

fa 
< 

Eh 
H 
U 


fa 

<: 
z 
o 

H 

CO 

fa 

S 

H 
Q 
I 

Z 

fa 

o 

fa 

HI 
> 

fa 
> 

o 


CM 


fa 
fa 

CJ 

H 
fa 


12 


Fig.  1.3. a.   An  example  of  two-dimensional  filtering 
(An  original  photograph)  [2] 
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Fig.  1.3.b. 


An  example  of  two-dimensional  filtering 
(Low  pass  filtered  version  of  Fig.  1.3a) 
[2] 
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Fig.  1 . 3c. 


An  example  of  two-dimensional  filtering 
(High  pass  filtered  version  of  Fig.  1.3a) 
[2] 
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nonrecursive  one  in  number  of  computations  and  storage 
requirements.   Hall  demonstrates  that  the  recursive  filter 
algorithm  constitutes  the  best  method  for  large  data,  i.e., 
it  is  the  fastest  and  cheapest. 

When  processing  two-dimensional  data  by  recursive  filters 
a  fundamental  problem  exists  due  to  inherent  feedback,  namely, 
the  problem  of  numerical  stability.   Since,  in  general,  two 
variable  polynomials  cannot  be  factored  into  a  product  of 
first  and  second  order  real  coefficient  polynomials  in  each 
of  the  variables,  it  is  difficult  to  solve  the  stability 
problem  for  two-dimensional  recursive  filtering.   Conse- 
quently, the  majority  of  papers  published  in  two-dimensional 
digital  filtering  deal  with  the  design  of  nonrecursive 
filters  which  are  inherently  stable  [6-10]. 

There  are  several  papers  discussing  two-dimensional 
recursive  digital  filters.   For  example  [11]  and  [12] 
formulate  two-dimensional  recursive  filters  by  the  z- 
transform  and  linear  difference  equations,  and  although 
they  investigate  problem  areas  related  to  stability  and 
realization,  the  majority  of  problems  remain  to  be  solved. 

The  most  important  applications  of  two-dimensional 
digital  filtering  of  digital  data  are  in  picture  processing 
and  geophysical  data  analysis.   Picture  processing  can  be 
categorized  into: 

a)  Digital  image  restoration  and  enhancement, 

b)  Computer  pictorial  pattern  recognition. 
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The  methods  of  image  restoration  and  enhancement  are 
applied  to  invert  degradations,  such  as  aberrations, 
atmospheric  effect,  scanning,  motion,  and  to  manipulate 
images  to  improve  viewing  phenomena  experienced  by  the  human 
eye.   Image  restoration  and  enhancement  has  been  used  with 
great  success  in  biomedicine,  i.e.,  in  extraction  of  quan- 
titative information  from  x-ray  films,  chromosome  counting, 
measuring  the  extent  of  arteriosclerosis  from  arteriograms, 
in  forensic  sciences  application,  i.e.,  fingerprint  image 
enhancement  for  automatic  classification,  and  in  astronomy, 
i.e.,  removal  of  turbulence  from  astronomical  photography 
[13]. 

Pictoral  pattern  recognition  by  digital  computer  has 
gained  great  importance  in  satellite  surveillance  and  mili- 
tary reconnaissance.   The  application  of  two-dimensional 
digital  filters  permits  the  separation  of  different  horizontal 
scans  on  magnetic  and  gravity  maps  in  geophysics  and  can 
be  used  equally  well  for  structural  and  topographic  maps 
or  for  any  other  type  of  factor  which  is  available  in  the 
format  of  a  planar  grid  [14]  .   Although  digital  filtering 
in  several  dimensions  is  gaining  importance  in  medicine, 
i.e.,  heart  volume  measurements,  in  fire  control  problems, 
and  to  analyze  complex  electronic  circuitry,  there  exists 
no  general  theory  concerning  structures,  analysis  and 
design. 

These  considerations  can  be  summarized  as  follows. 
Since  one-dimensional  recursive  digital  analysis  and  design 
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techniques  depend  strongly  on  the  fundamental  theorem  of 
algebra  and  on  extensive  utilization  of  continuous  design 
theory,  it  represents  a  special,  i.e.,  a  non-generalizable 
field  in  the  area  of  N-dimensional  recursive  digital  filtering. 
Multi-dimensional  digital  filtering  is  well  developed  in  two- 
dimensions  but  is  based,  due  to  the  absence  of  recursive  theory 
analysis  and  design  tools,  predominantly  on  nonrecursive 
filtering  schemes.   Because  of  their  inherent  advantages,  much 
of  the  current  research  is  directed  towards  recursive  filtering 
algorithms  and  techniques. 

B.   AREA  OF  INVESTIGATION  AND  OBJECTIVES 

One  area  of  investigation  of  this  thesis  is  to  see  how  the 
semi-infinite  straight  line  approximation  technique  of  one- 
dimensional  filters  (based  upon  the  log  modulus/log  frequency 
or  decibel  vs.  log- frequency  plots)  can  be  extended  to  two 
(or  higher)  dimensional  recursive  filters.   The  advantages  of 
log-modulus  approach  is  that  a  transfer  function  in  factored 
form  becomes  a  linear  sum  of  terms  when  the  logarithm  of  its 
magnitude  is  taken,  and  each  of  these  terms  can  be  approxi- 
mated by  semi-infinite  straight  line  segments  in  the  one- 
dimensional  case,  and  semi-infinite  planes  in  the  multi- 
dimensional case,  when  expressed  in  terms  of  log-frequency. 

The  following  are  specific  objectives: 

1.   To  postulate  typical  (canonic)  factors  for 
two-dimensional  filters  and  to  investigate  their  properties 
in  the  log-modulus/log  frequency  domain. 


18 


2.  To  develop  an  approximation  technique  for  the 
design  of  two-dimensional  recursive  filters  to  meet  given 
specifications  in  the  frequency  domain. 

The  frequency  response  of  two  dimensional  filters  is 
inherently  a  volume  in  the  space  of  magnitude,  log  cj,  and 
log  CO2  •   This  volume  is  the  given  or  desired  specification 
of  the  filter,  which  in  the  separable  case  can  be  approxi- 
mated by  the  combination  of  semi-infinite  planes. 

In  this  thesis  we  investigate  how  we  can  approximate 
this  volume  by  semi-infinite  planes  as  an  extension  of 
the  approximation  used  in  one  dimension,  i.e.,  by  semi- 
infinite  lines.   After  determining  the  planar  equations 
which  perform  the  given  specification  of  the  filter,  we 
investigate   ways  to  get  a  transfer  function;  first,  in  the 
continuous  frequency  domain  and  second  for  two-dimensional 
recursive  equations  in  the  discrete  variable  domain. 

The  results  of  this  study  will  classify  those  charac- 
teristics that  can  be  achieved  by  cascading  simple  canonic 
factors,  and  indicate  limitations  that  are  imposed  by 
trying  to  achieve  a  particular  frequency  characteristic 
(such  as  the  fan  filter)  in  as  simple  form  as  possible. 

The  approach  taken  here  is  to  investigate  the  problem 
in  the  continuous  domain  first,  and  then  study  the  trans- 
formation from  continuous  frequency  (differential  equation) 
domain  to  discrete  time  (sampled  data)  domain. 
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C.   PREVIEW  OF  RESULTS 

We  have  proposed  a  geometric  method  to  design  two- 
dimensional  recursive  digital  filters.   This  is  a  new 
method  which  is  relatively  simple  and  easy  with  respect  to 
other  design  methods.   We  have  investigated  properties  of 
some  typical  (canonic)  factors  and  considered  their  stability 
properties.   In  Chapter  II  we  go  through  the  basic  properties 
of  two-dimensional  filtering.   In  Chapter  III  we  provide 
background  in  terms  of  planar  geometry  for  the  discussion 
in  subsequent  chapters.   In  Chapter  IV  we  investigate  the 
properties  of  typical  factors  and  their  stability  in  continu- 
ous variable  domain.   The  new  design  technique,  including 
stability  considerations,  starting  with  the  continuous 
frequency  domain  and  extending  to  two-dimensional  recursive 
filters  in  the  discrete  domain,  is  presented  in  Chapter  V. 
The  results  are  summarized  and  some  comments  for  further 
study  are  made  in  Chapter  VI.   It  should  be  noted  that  the 
results  presented  apply  to  quarter  plane  filters,  as  defined 
in  Chapter  V. 
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II.   FUNDAMENTALS  OF  TWO-DIMENSIONAL  RECURSIVE  FILTERING 

A.   INTRODUCTION 

There  are  many  signals  that  are  inherently  two-dimensional 
for  which  two-dimensional  signal  processing  techniques  are 
required.   The  most  common  examples  are,  of  course,  images 
in  which  the  two  variables  are  the  spatial  coordinates. 
Image  processing  [15],  [16]  plays  an  important  role  in  many 
areas  of  scientific  and  technical  research.   Some  examples 
are  satellite  imagery  radiography,  radar,  and  biomedical 
images  such  as  acoustical  holograms,  medical  x-rays,  and 
electron  micrographs  [2].   Not  all  two-dimensional  data 
come  from  an  image.   In  prospecting  for  gas  and  oil,  a 
common  procedure  is  to  monitor  the  seismic  signals  produced 
by  detonating  an  explosive  charge  [17] .   An  array  of  geo- 
phones  situated  along  a  line  radiating  in  a  vertical  or 
horizontal  direction  from  the  point  of  explosion  is  used 
to  record  these  signals.   In  this  case,  the  two  variables 
are  distance  (from  the  point  of  explosion  measured  along 
the  line)  and  time  (the  time  at  which  the  seismic  signal 
reaches  a  given  geophone  or,  equivalently ,  a  given  distance 
from  the  explosion) . 

Although  two-dimensional  signals  may  be  processed  by 
one-dimensional  systems,  it  is  generally  preferable  to  con- 
sider using  two-dimensional  systems.   Many  of  the  basic 
ideas  of  one-dimensional  signal  processing  may  be  readily 
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extended  to  the  two-dimensional  case.   There  are  some 
very  important  concepts  of  one-dimensional  systems,  however, 
that  are  not  directly  extentable  to  two  dimensions.   It 
is  the  goal  of  this  chapter  to  discuss  the  basic  ideas  and 
techniques  of  two-dimensional  signal  processing  and  to 
illustrate  them  into  the  context  of  two-dimensional  filter 
design.   A  more  detailed  discussion  of  this  theory  plus 
an  overview  of  recent  work  in  two-dimensional  filtering  is 
contained  in  the  excellent  survey  paper  by  Mersereau  and 
Dudgeon  [18] ,  and  H.  Chang  and  J.K.  Aggarwal  [19] .   Other 
useful  references  are  the  first  few  chapters  of  [1]  and 
chapter  7  of  [2] . 

B.   DEFINITIONS 

1.   Two-Dimensional  Signals 

One-dimensional  signals  are  functions  of  a  single 
variable  and  two-dimensional  signals  are  functions  of  two 
integer  variables.  Consider  the  two-dimensional  sequence 
x(m,n)  where  m  and  n  are  integer  variables,  shown  in  figure 
2.1.  As  in  the  one-dimensional  case,  the  notation  x(m,n) 
is  often  short  hand  for  a  sampled  version  of  a  continuous 
two-dimensional  signal  x(s,t):   i.e., 


x(m,n)   =   x(mT,,nT2)   =   x(s,t) 


(2.D 

s=mT1,t=nT2 


One  interpretation  of  x  is  a  function  which  assigns  a 
(generally  complex)  number  to  each  integer  ordered  pair 
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-7y-  it/aM 


* 


Fig.  2.1.   Pictorial  representation  for  two- 
dimensional  discrete  array 


X  P £$ flit  \-$l> J 
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(m,n) .   The  function  x  is  not  defined  when  either  or  both 

of  its  arguments  (m,n)  is  not  an  integer.   Such  a  signal 

will  be  interchangeably  referred  to  as  either  a  two-dimensional 

array,  a  two-dimensional  sequence  or  simply  as  a  two-dimensional 

signal . 

Some  useful  two-dimensional  sequences  are  defined 
below  and  shown  in  figures  2.2  to  2.4.   These  include: 

1.  Digital  Impulse  or  Unit  Sample 

1    m  =  n  =  0 
u  (m,n)   =  \  (2.2) 

0     elsewhere 

2.  Digital  Step 

( 

I  1    m,n  >_  0 
u  ,(m,n)   =  \  (2.3) 

0     elsewhere 

3.  Exponential 

m  n        .  n 
a,  a~    m,n    0 

x(m,n)   =      X  (2.4) 

0        otherwise 

4.  Sinusoid  (Complex) 

j  ((jo^m+oj^n) 
x(m,n)   =   e  -°°  <  m,n  <  +=°  (2.5) 

As  seen  above,  the  two-dimensional  step  is  related  to  the 
two-dimensional  impulse  by 
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■+ 


+ 


+-H-4 


-t 


■f-  f-  -h-h 

■/  -h-h  -h  +-  /-■/-  +-/-++- 
-h-  -h  ■/-  f-  -H  f  4    A  -h   A  4  +■ 

f  +■  4  +  f  4  4  4  -h  +  + 
4-  4-  4-  -hif-  4  4-  4  -h-f-  -h 


-4  -h  -i-4-  f-  -f-  +  -h-h-h  4  i- 
-h -h  4- -h  +-  -f-f-  4  1-+-  j-4- 

+  j-4  4-  f  -h-  4-  ■/-  4  -h  -h 
i--l-4--/-4-  -/-4-4--h+- 
-h4—h   4-/-  4-  4  +■  4  4 


I     /  / >  rn. 


Fig.  2.2.   Digital  Impulse 
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-     I  *  *  -f  4 

f   +■  +  -t-_zf-  X  XXX  X 

-h-h  +  +  f  XXXXXX 
-h  +  +  4-4- XX  xxxx 
-h-h-h  h-h.XX  XXXX, 
-h-t-  4- 4- -4  .XXX  XXXX 

+.+  +  +  i.XiXXxxx 

/  /  /  /  /  /  y  )  }  y  }  )  >" 

4--h4-4-4-4-4-4-i--h4-  4- 

-h  -h  4-  4-  f  -4  4   h^4   4 

4--A  4-  4  4  4  4  4  +  4--h 
-h  4  -h-j-h-  4-  4-  -h  4-  -h 
+  +  -h  -/  4  -L-  -h  4  4  4 


Fig.    2.3.       Digital    Step 
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-t-  -h  H-  -h  -j-  -J-  -f-    ■f-  f-   -f- 

-*-  +  +  ++•  f--*-    ff  f 
-/-  -f-  ■+   4-+   f  4-    4-f-   *■  *■ 

-f-   +  +  -h  -h  f 1    *   f-   *-* 

■t-  -t  +  -h  HfA-  J-  j-  -f-  f-  + 


m 


4  4  -h  -h  4  4  4-  +  -f-  -h  -h  +■ 
4-  +  4  +  +  4- 4- +  -h  +- -h  -t- 

+-4-4  -h  4  4-  4  4-  4  4  f- 
-4  4-  7-/  -h  -h  A  -h  +~  4- 


Fig.    2.4.       Exponential 
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m     n 
u_1(m,n)   =11   uo(m1,n1)        (2.5a) 


m,  =-<»  m,  =-oo 


2 .   Two- Dimensional  Systems 

A  two-dimensional  discrete  system  is  defined  as 
the  unique  transformation  T  that  maps  an  input  sequence 
(x(m,n)}  and  the  initial  condition  sequence  {s(m,n)}  into 
an  output  sequence  {y(m,n)} 

y(m,n)   =   T[{x(m,n) },{s(m,n) }]  (2.6) 

In  signal  processing  applications,  we  assume  zero  initial 
conditions.  With  this  assumption,  the  above  relation  may 
be  shortened  to 

y(m,n)   =   T[{x(m,n)}]  (2.7) 

The  system  characterized  by  T  is 

(a)  linear  if 

T[{ax1(m,n)  + bx2 (m,n) }]   =   aT [ {x1 (m,n) } ]  + bT [ {x2 (m,n) }] 

(2.8) 
for  arbitrary  constants  a  and  b. 

(b)  shift  invariant  if 

y(m-m1,n-n1)   =   T  [  {x  (m-ir^ ,  n-n1)  }  ]  (2.9) 

for  all  m,  and  n, . 

28 


A  hardware  and  software  system  which  is  equivalent  to  the 
transformation  T  of  (2.6)  is  called  a  digital  filter.   The 
class  of  complex-valued  sequences  which  are  square  summable 
will  play  the  role  of  input  and  output  signals  of  the 
digital  filter.   Linearity  and  shif t-invariance  are  indepen- 
dent properties  of  a  system  in  that  neither  property  implies 
the  other.   Thus  the  system 

T[x(m,n)]   =   h(m,n)  x(m,n) 

is  linear  but  not  shift  invariant  and  the  system 


2 

T[x(m,n)]   =   x  (m,n) 


is  shift  invariant  but  not  linear 

For  two-dimensional  discrete  linear  systems,  the 
impulse  response  is  defined  as  the  response  of  the  system 
at  (m,n)  to  a  two-dimensional  digital  impulse  input  at 
(m, /n,)  with  zero  initial  conditions,  i.e., 

h(m,n;m,,n, )   =   T [ {uq (m-m1,n-n1) } ]      (2.10) 

where  u  (m,n)  denotes  the  digital  impulse.   In  particular, 
the  impulse  response  of  a  two-dimensional  discrete  LSI 
(Linear  Shift  Invariant)  system  is 
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h (m,n;m, ,n, )   =   h (m-m, ,n-n, ) 


or 


h(m,n)   =   T[{uo(m,n)  }]  (2.11) 


The  response  of  a  two-dimensional  LSI* system  with  zero 
initial  conditions  is  completely  characterized  by  its 
impulse  response.   In  fact  using  the  identity 


x(m,n)       =  I  x(m, ,n ,  ) u    (m-m, ,n-n1 ) 


l'"l'"o% V      "1' 


m1,n1=-« 


I  x(m-m1,n-n1)u    (m1,n1) 


m^n^-00 


we   get 


y(m,n)      =      T[{  x (m, ,n1 ) u    (m-m, ,n-n, ) } ] 


rarnr 


l  x(m1,n1)T[{u0(m-m1,n-n1) }] 


m1/n1=-c 


x (m, ,n, ) h (m-m, ,n-n^) 


ml'nl=~°°  ' 


*LSI  =   Linear  Shift  Invariance 
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y(m,n)   =     £     x  (m-m^n-r^)  h  (m,  ,n,) 


m1/n1=-< 


=   x(m,n)  *  h(m,n)  (2.12) 

where  "*"  denotes  two-dimensional  convolution.   In  other  words, 
for  LSI  systems  the  basic  convolution  theorem  is  valid. 
3.   Causality,  Separability,  Stability 

In  discussing  various  features  of  two-dimensional 
discrete  systems,  it  is  convenient  to  have  spectral  nota- 
tion [19]  to  represent  certain  regions  of  the  spatial  domain. 

1.   Let  S(a,3)  denote  a  sector  with  the  angular 
interval  (a, 3).   In  polar  form  S(a,3)  may  be  described  by 

S(a,3)   =   ( (r,6) |r  > 0,a  <  9  <  3)  (2.13) 

where  (a, 3)  are  measured  in  radians.   In  particular,  we 
give  the  following  notation  for  the  sectors  frequently 
encountered: 

S++      =      S[0,tt/2];         S-+      =      S[Tr/2,ir] 

S—      =      S[-tt,-tt/2]  ;       S+-      =      S[-tt/2,0]        (2.14) 

S*+   =   S[0,tt)  ;      S*-   =   S[-tt,0) 
where  "["  and  "]"  imply  including  the  boundary  and  ")"  implies 


X 


not  including  the  boundary. 
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Notice  that  S++,  S-+,  S — ,  and  S+-  represent  first,  second, 
third,  and  fourth  quadrants,  respectively,  and  S*+  and  S*- 
represent  half-planes  which  are  symmetric  to  each  other  with 
respect  to  the  origin.  (See  fig.  2.5.) 

2.   For  a  two-dimensional  sequence  (x(m,n)},  the 
set 


Ev      =      { (m,n) |x(m,n)  ?   0}  (2.15) 


is  called  the  support  of  (x(m,n)}. 

3.  Let  {h(m,n)}  be  the  impulse  response  of  a  two- 
dimensional  discrete  LSI  system,  and  let  E,  be  the  support 
of  (h(m,n) }. 

Then,  a  two-dimensional  LSI  system  is  said  to  be 
causal  if  E,  is  the  set  S++,  and  semi-causal  if  E,  is  the 
set  S*  +  .   They  are  shown  in  Fig.  2.5. 

A  two-dimensional  discrete  LSI  system  is  said 
to  be  separable  if  its  impulse  response  can  be  factored  into 
a  product  of  one-dimensional  responses  (Fig.  2.6);  i.e., 

h(m,n)   =   hx(m)  h2  (n)  (2.16) 

If  the  equation  (2.16)  is  not  satisfied,  the  filter  is  said 
to  be  nonseparable.   The  advantages  of  separable  filters  is 
that  two  dimensional  convolution  (2.12)  may  be  carried  out 
as  a  sequence  of  one-dimensional  convolutions.   This  can  be 
seen  by  rewriting  (2.16)  as  follows. 
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fcm 


Fig.  2.5.   Support  Area  of  Impulse  Response 

(a)  Causal 

(b)  Semi-causal  S. 

(c)  Semi-causal  S*_ 
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X(«r«.) 

H,(-0 

H2(«.) 

) 


1'  2 


H(z)   =   H1(z1)H2(z2) 


Fig.  2.6.   Separable  Digital  Systems 
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y(m,n)   =      j»     h1  (m^)  h2  (r^)  x  (m-m-^n-n-^        (2.17a) 


m,  ,n  =-<» 


I        h1(m1) [  I        h2(n1)x(m-m1/n-n1) ]  (2.17b) 


m  =-oo        n,=-oo 


I        h1(m1)  a(m-m.,n)  (2.17c) 


ml=_c 


where  a  (m-m,  ,  n)  is  a  sequence  of  one-dimensional  convolutions 
(obtained  by  evaluating  the  terms  inside  the  bracket  of 
(2.17b)  for  each  fixed  value  of  m, ) .   Equation  (2.17c)  shows 
that  y(m,n)  may  be  obtained  by  a  second  sequence  one-dimensional 
convolution. 

If  both  the  input  sequence  x(m,n)  and  the  filter 
impulse  response  h(m,n)  are  separable,  then  it  is  readily 
seen  that  the  output  sequence  y(m,n)  is  also  separable. 
In  this  case  we  obtain  the  result  [using  (2.12)  and  (2.16)] 

00  00 

y(m,n)   =   [  I      h±  Cm^b  (m-m1)  J  I  I      h2  (n1)  c  (n-i^)  ]  (2.18a) 
m  =-«  n  =-oo 

=   a(m)3(n)  (2.18b) 


where 


x(m,n)   =   b(m)  c(n) 
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A  two-dimensional  LSI  system  is  said  to  be  (BIBO) 
stable  if  every  bounded  input  sequence  produces  a  bounded 
output  sequence.   This  condition  is  satisfied  if  and  only 
if  the  two-dimensional  impulse  response  of  the  system  is 
absolutely  summable,  i.e., 

00 

I  |h(m,n)  |   <   c°  (2.19) 

m,  n=-°° 

As  in  the  one-dimensional  case,  (2.19)  can  be  shown  to  be  a  <e- 
necessary  and  sufficient  condition  for  stability  [20].   One 
problem  with  (2.19)  is  that  it  can  be  quite  difficult  to 
evaluate  for  an  arbitrary  h(m,n). 

4 .   Two-dimensional  Difference  Equations 

As  in  the  one-dimensional  case,  two-dimensional  LSI 
systems  can  often  be  described  by  a  two-dimensional  linear 
constant-coefficient  difference  equation  relating  the  output 
y(m,n)  of  the  filter  to  its  input  x(m,n).   The  most  general 
form  for  such  a  difference  equation  is  shown  by 


I  a(k1,k2)x(m-k1,n-k2)  -     I  b  (  ^^  y  (m-Jl^n-j^)  =  0 

(k1,k2)eRa  (i1,l2)eRh 

(2.20) 


wh 


ere  R  and  R  are  finite  sets  of  spatial  grid  points, 


called  the  input  and  output  mask,  respectively,  and   {a(k1,k2)} 
and  (bUwJU)}  are  the  set  of  constant  coefficients  that 
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characterize  the  particular  filter.   Generally,  a  difference 
equation  has  a  family  of  solutions  as  is  the  case  with  a 
differential  equation.   The  difference  equation  is  said  to 
be  recursive  if  there  exists  a  sequence  of  computations  which 
yields  the  output  sequence  serially  from  the  input  sequence 
and  that  portion  of  the  output  sequence  which  is  already 
computed/  including  the  initial  condition.   An  implicit 
ordering  of  grid  points  in  the  spatial  domain  is  assumed  for 
the  serial  computation  of  the  output  sequence. 
There  are  two  cases  of  (2.20): 

1.  R  ,  K     in  the  S++  (first  quadrant),  b(0,0)  ^  0 
and  y(m,n)  =  0  for  outside  of  S++  (zero  initial  conditions), 
then  (2.20)  represents  a  causal  system. 

2.  If  R  ,  R.  in  S*+,  b(0,0)  ?   0  and  y(m,n)  =  0  for 

a   o 

outside  of  S*+  (zero  initial  conditions),  then  (2.20) 
represents  a  semi-causal  system. 

In  both  cases  we  may  assume  b(0,0)  =  1  without  loss  of 
generality  and  write  (2.20)  as 


y(m,n)   =      £     a  (k-^k^  x  (m-k-^n-k^ 
(krk2)£Ra 


I  b(l1,!L2)Y(m-l1,n-l2) 

(2.21) 


(llfl2) eRb-(0,0) 


It  is  clear  that  (2.21)  is  a  recursive  equation,  and  the 
implicit  ordering,  for  example  may  be  given  as 
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{(0,0); (1,0), (0,1),  (1,1); (2,0) ,(2,1) , (0,2) ,(1,2); (2,2);...} 

for  a  causal  system,  and 

{(0,0), (0,1) , (0,2) ,...; (-1,1) , (0,1) , (1,1) ;...;... 

...  (-1,2)  ,  (0,2)  ,  (1,2)  ,..  .;  .  ..} 

for  a  semi-causal  system. 

If  {b(£1,£2)}  =  {uQ(£1,£2)}  the  equation  (2.21) 
becomes 

y(m,n)   =      I  a(k1,k2)x(m-k1,n-k2)    (2.22) 

(k1,k2)£Ra 

so  that  the  output  sequence  is  a  weighted  moving  average  of 
the  input  sequence.   The  impulse  response  of  (2.22)  is 
unique  and  is  obviously  {a(k,,k2)}.   Such  filters  are  called 
finite  extent  impulse  filters  (FIR)  because  the  impulse 
response  has  only  finitely  many  nonzero  terms.   FIR  filters 
obviously  are  always  BIBO  stable  because   a(k,,k2)   is  always 
absolutely  summable.   (It  has  only  finitely  many  nonzero 
terms . ) 

If,  on  the  other  hand,  the  sequence  (b(£,,£2)} 
has  several  nonzero  terms,  then  the  impulse  response  will 
generally  have  infinitely  many  terms.   We  will  say  that  such 
filters  are  infinite-extent  impulse  response  filters  (IIR)  . 
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Such  filters  may  be  BIBO  unstable,  and  so  it  is  necessary 
to  devise  a  stability  test  for  such  filters. 

FIR  filters  have  the  advantages  that  they  are 
easier  to  design,  that  controlling  the  phase  response  is 
easier  (in  particular,  it  is  easy  to  produce  filters  which 
have  a  linear  phase  response) ,  and  that  since  the  input- 
output  relationship  is  a  finite-extent  convolution,  Fast 
Fourier  Transform  (FFT)  techniques  can  be  used  to  speed  the 
computation  of  the  output  sequence.   On  the  other  hand, 
recursive  (IIR)  filters  generally  require  less  high  speed 
storage  and  usually  require  less  computation  time  than  FIR 
filters.   Comparison  of  FIR  and  IIR  filters  is   contained 
in  [5]. 

5 .   Two-dimensional  z-tranform 

A  useful  mathematical  tool  for  representation  of  a 
sequence  x(m,n)  is  its  two-dimensional  z-transform.   Given 
a  sequence  x(m,n),  the  z-transform  of  this  sequence  is 
defined  to  be: 

< 

CO        oo 

Z[{x(m,n)}]   4  X(z1,z2)   =     I  £  x  (m,n)  z^   (2.23) 

m=-°°  n=-°° 

where  z,  and  z2  are  independent  complex  variables.   The  most 
useful  property  of  the  z-transform  is  the  fact  that  the  z- 
transform  of  the  convolution  of  two  sequences  is  the  product 
of  their  z-transforms ,  i.e.: 
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Z[{x(m,n) }*{h(m,n) }]   =   X(z1,z2)H(z1,z2)  (2.24) 


Therefore,  the  input-output  relationship  of  an  LSI  filter 
may  be  expressed  in  the  z-transform  domain  as: 

Y(z1,z2)   =   H(z1,z2)X(z1,z2)  (2.25) 

where  H(z-,,z2),  the  z-transform  of  the  impulse  response 
is  called  the  transfer  function. 

A  two-dimensional  transfer  function  evaluated  on 
the  unit  circle  is  called  the  frequency  response  of  the 
system,  i.e., 

jw,       jco9  ju>,       jo>9  jui,       joo9 

Y(e         ,e      *)      =      H(e      x,e        )X(e      x,e      z)  (2.26) 

when  a  two-dimensional  LSI  system  is  describable  by  the 
difference  equation  of  (2.20),  and  we  take  a  z-transform    >' 
of  this,  it  follows  that 


-k   -k2 
[      I  a(k1,k2)z1   z2   ]X(.z1,z2) 


(kx,k2) eRa 


-  [     I  b(£1,£2)z1   z2  Z]Y(zL,z2)   =   0 

(ilfl2)  el^ 
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or 


H(z1,z2)   =   Y(z1,z2)/X(z1,z2) 


'kl  ~k2 
I  a(k1,k2)z1   z2 


(kl'k2)eRa 


-A,  -l0 
I  b(£1,£2)z1  xz2  z 

( *■  -i  /  & ^ )  s^j-. 


=   A(z1,z2)/B(z1,z2)  (2.27) 

where  A  and  B  are  the  z-transf orms  of  (a(k,,k2)}  and 
{b  (£■•*  &2)  }  respectively.   In  particular  when  (2.27)  repre- 
sents a  causal  or  semi-causal  system,  i.e.,  (2.27)  is  the 
z-transform  of  (2.21),  we  can  give  explicit  descriptions  for 

R  and  R,  in  each  case.   For  the  causal  case  we  define 
a      d 

(Fig.  2.5a) 

Ra      =      *  (ki'k2>  I  °  ikl  iMl'    °ik2-M2} 

(2.28) 

Rj^      =      (  (llfi2)  |  0  <_l1  <_Nlf    0£^2-N2^ 
For   the    semi-causal   case   we   define    (Fig.    2.5b) 

Ra      =      {  (ki'k2}  1°  ikl-rVk2  =  °;"Mct  ikliM6'°  <  k2  -M2  } 

(2.29) 

Rb      =       {(£1,£2)|0<^llN3,£2  =  0;-Nal£1<Na,Q<£2lN2} 
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and  in  both  cases  b(0,0)  =  1  is  assumed.   The  inverse 
z-transform  is  defined  as 


x(m,n)       =      j  j        /      X(z,,z    )  z™   1    z^        dz,dz„     (2.30) 

(2njr    n        „  12      12  12 

Cl      c2 

where  c^  and  c2  are  suitable  closed  contours  in  the  z,  and 
z2  planes. 

The  two-dimensional  Fourier  transform  of  a  sequence 
is  defined  as 


DO).,  3u? 

X(e    ,e  *)      =  F[x(m,n)] 


X(z1,z2) 


2 
(z1,z2) eF 


-  j  ( co-j  m+cu^n) 
£    x(m,n)e     x    *  (2.32) 

mr  n=-°° 

For  any  given  sequence  the  set  of  (z,,z2)  for  which  the 
z-transform  converges  is  called  the  region  of  convergence. 
Uniform  convergence  of  the  z-transform  requires  that  the 
sequence  be  absolutely  summable,  i.e. 

oo 

I         |x(m,n)  |  iz1|"m|z2|"n   <   =o        (2.32a) 
m,  n=-°° 
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Similarly,  the  Fourier  transform  X(e    ,e    )  converges 
uniformly  to  a  continuous  function  of  co,  and  w_  if 


I  |x(m,n)  |   <   oo  .  (2.32) 

m,n=-°° 


A  power  series  of  the  form  of  (2.23)  is  known  as  a  two- 
dimensional  Laurent  series  in  the  theory  of  complex  functions 
in  two  variables . 

6 .   Two-dimensional  DFT 

Let  {x(m,n)}   be  a  finite  support  sequence  of  the 
size  MxN.   Then  the  two-dimensional  discrete  Fourier 
transform  (DFT)  of  (x(m,n)}   is  defined  as 


M~1N~X  -^(Rki  +  £k2> 

X(k1#k2)   =11    x(m,n)e     M  L       N  z     (2.33) 

m=0  n=0 


for  0  <_  k.  <  M-l  and  0  <_  kp  <  N-l.   Similarly,  the  inverse 
discrete  Fourier  transform  (IDFT)  is  given  by 
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X(m,n) 


M-l 

-  I 


MN 


k,=0 


N-l 
I 

k2=o 


x(klfk2)e 


k,    k2 

+:2TT(irm  +  irn) 


(2.34) 


Note  that  the  DFT  corresponds  to  sampling  the  Fourier  trans- 
form of  MxN  points,  i.e., 


X(kx,k2) 


=   x  (e 


Ua>i   300. 
/e 


kl       k2 

a)l=27TlT'a)2=27T-N- 


(2.35) 


]W. 


30). 


More  generally,  if  X(e    ,e    )  is  the  Fourier  transform 
of  (x(m,n)}   whose  size  is  greater  than  MxN,  then  the 
IDFT  of  {X(k1 ,k  )}  will  result  in  aliasing.   In  other  words, 
the  IDFT  of  (X(k,,k2)}  will  yield  one  period  of  (x(m,n)} 
where 


x(m,n)   = 


m1,n1=-<» 


x  (m+m,M,n+n,N) 


(2.36) 
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C.   STRUCTURES  OF  TWO-DIMENSIONAL  FILTERS 

1.   One-dimensional  Digital  Filter  Realization 

One-dimensional  digital  filters,  in  general,  can 
be  shown  in  the  z-transform  domain  as  follows: 

y 


1=0  Y(zi) 

H<zl>  =  ~ =  xuJt         (2'37) 

1  -  I      am  Zl 
m=l 


or ,  in  the  time  domain 


M  L 

y(n)       I      a   y(n-m)  +   J   b  x(n-l)  (2.38) 


m 
m=l  1=0 


The  realization  implementing  equation  (2.38)  directly,  is 
known  as  the  direct  form  1  and  shown  in  Figure  2.7. 

Equation  (2.37)  can  be  rewritten  in  a  slightly 
different  form  by  introducing  a  new  variable,  W(z1) ,  such 
that 


W(z,)    Y(z,) 

H/z  )   =   L_  .   .  xs  (2.39) 

H(Z1J      X(zx)    W  z1 


The  corresponding  set  of  difference  equations  consists  of: 
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M 
w(n)   =   x(.n)  +  I      am  w(n-m)  (2.40) 

m=l 


and 

L 
y(n)       I      b£  w(n-£)  (2.41) 

1=0 

Equations  (2.40)  and  (2.41)  are  realized  in  the 
direct  form  2  as  shown  in  Figure  2.8  which  can  be  visualized 
as  a  cascade  realization  of  two  digital  filters,  realizing 
the  denominator  and  numerator  polynomials,  respectively. 
Figure  2.8  can  be  redrawn  as  shown  in  Figure  2.9.   The 
resulting  realization  is  called  direct  form  3  or  canonic 
form,  since  it  represents  a  structure  with  the  number  of 
delays  equal  to  the  order. 

At  this  point,  an  additional  realization  will  be 
introduced,  the  direct  form  4,  which  combines  the  charac- 
teristics of  the  direct  forms  1  and  2  having  only  two 
entries  in  each  summer,  which  corresponds  directly  to 
hardware  implementation  and  by  realizing  the  numerator  and 
denominator  polynomials  separately  in  cascaded  form 
(Figure  2. 10) . 

There  are  several  other  methods  to  realize  a  one- 
dimensional  M   order  digital  filter.   For  example,  the 
cascade  form  of  first  and  second  order  sections  H(z,),  as 
shown  in  Figure  2.11,  where 
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H(zX) 


2         b£    Zl 
£=0 


-£ 


M 


1    "      I      a      z. 

L        ml 

m=l 


-m 


X(n) 


W(n) 


I 


«— 2 

H 

«,  r 

H 

rf                   1 
^           3              V 

7-1 

Fig.    2.8.       Direct   Form    2 
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H(z1) 


£=0         *       L 


M 


i  -     y    a    z" 

L        ml 
m=0 


-m 


Y(-) 


Fig.    2.9.       Direct   Form    3 
(Canonic) 
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H(z1)   =   aQ   n   Hi(z1),   *  =  ^        (2.42) 
i=l 


Each  component  H^(z^)  can  be  realized  in  one  of  the  above 
outlined  forms. 

Partial  fraction  expansion  methods  applied  to 
equation  (2.37)  lead  to  a  parallel  form  realization  of 
first  and  second  order  sections  (Figure  2.12),  i.e., 

k 
H(z1)   -   c+  J      Hi(z1),   k=  [^+1]  (2.43) 

i=l 

Other  forms  of  realization  include  hybrid  structures, 
i.e.,  parallel-cascade  forms,  the  transpose  configuration, 
which  can  be  obtained  for  all  previously  outlined  structures 
by  reversing  the  direction  of  signal  flow  and  by  interchanging 
all  branch  and  summing  modes,  wave-digital  and  continued 
fraction  expansion  filters. 

The  direct  forms  discussed  above  and  the  series/parallel 
arrangement  of  lower  order  sections  are  by  far  the  most 
often  used  ones  in  computer  simulation  and  in  digital 
hardware. 

2.   Two-dimensional  Digital  Filter  Realization 

Two-dimensional  Quarter  plane  digital  filters,  in 
general  form,  can  be  shown  as 
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HW  =  C+2H}U,) 


Fig.  2.12.   Parallel  Form 
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£  =0  £  =0    X  Z 
H(z,,z  )   =   — - (2.44) 

Ml    M2 

-m..   -m 

1  ~   I     I    ammz,    z0 
L     *<     m.iru   1     2 

ml=0  m2=° 

The  numerator  of  the  two-dimensional  transfer  function  (2.44) 
can  be  written  as  a  one-variable  polynomial  in  z~,    weighted 
by  coefficients  which  are  functions  in  z,,  i.e., 


L2 

Ll 

I 

(  I 

«,,=0 

J..=0 

L2 

'SV 

£1   _il2 


z2     (2.45) 


"£2 
:„  (2.46) 


The  same  reasoning  applied  to  the  denominator  polynomials 
leads  to 


M     M 

-m,   -itu 

D(zrz2}    =    1  -     ^  (    l    ^2  Zl    )Z2 

m,=0  m.=0 

(2.47) 

M2 

=   1-1  (V.  <zl»»Z2m      (2-48) 
m2=0 
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To  obtain  the  direct  form  3  representation,  a  new  function, 


W(z1z2)  is  introduced: 


W(z,,z9)    Y(z,,z9) 

H(ziz2>    -   xi^Tzf)  ■  w(^fr  (2-49) 


where  each  factor  can  be  written  using  equation  (2.46)  as 


W(zx,z2)   = 


*2 


£2-0 


/NI1£2(Z1))Z2    *<V*2J 

(2.50) 


and  equation  (2.48)  as 


M2 


-m 


2 


Y(z1,z2)   =   W(z1,z2)  +    I    (DM  m  (z1))z2  Z  Y(z1z2) 

m2=0    ^  2 

(2.51) 

The  realization  of  equations  (2.50)  and  (2.51)  is  shown  in 
Figure  2.13  using  compact  notation,  and  Figure  2.14  where 

NL  l    (zl^  anc^  DM  m  ^zl^  are  ^mPlemented  f°r  each  (L-,^2) 
and  (nunu)  by  one-dimensional  direct  form  4  realizations. 

The  general  form  of  two-dimensional  causal  digital 
filter  (2.44)  can  be  rewritten  by  introducing  a  new  variable, 
W(z-,z2)  or  in  the  previous  function,  then, 

W  ( z ,  z  « )     Y  ( z ,  ,  z  9 ) 

H<VZ2>   "   X(Zl,z2)  •  W(z',z2)         (2"52' 
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Each  factor  can  be  written  as : 


Ll   L2 

-I  -I 

Y(ll'»2J   =     I  I      W  Zl  1    z2  2  WCzrz2) 

£1=0  £2=0 


(2.53) 


and 


M       M 

-m   -m« 
W(zrz2)   =-  XCVz2)  +   I  I      a     zx   z2   W(zrz2) 

m1=0   m2=0    -1- 


(m1+m2)^0)  (2<54) 


The  equations  (2.53)  and  (2.54)  can  be  written  as  a  differ- 
ence equation  as  follows, 


Ll     L2 
y(m,n)   =     I  I      bz    z      W(m-2.1,n-£2)      (2.55) 

ll=0      £2=0    1    2 


and 


Ml   M2 

w(m,n)   =  x(m,n)  +   I  I  am  m  w ^m"m2 ' n~m2 *      (2.56) 

m1=0  m2=0 


m1,m2  ^  0 
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Without  loss  of  generality,  L.  is  chosen  equal  to  M.  for 

all  i.   The  structure  realizing  equation  (2.55)  and  (2.56) 

is  shown  in  Figure  2.15.   It  is  noted  that  the  number  of 

delays  in  Figure  2.15  equals  the  order  of  H(z..,z  ).   The 

structure  is,  therefore,  by  definition,  canonic. 

The  direct  form  4  (Fig.  2.10)  realization  of  H(z.,z2) 

is  shown  in  Figures  2.16  and  2.17,  using  compact  notation  and 

one-dimensional  direct  form  4  (Fig.  2.10)  implementation  of 

each  N   «   and  D     ,  respectively. 
L1^2       Mlm2 

D.   STABILITY 

Two  major  problems  in  the  design  of  a  recursive  filter 
are  approximation  and  stability.   Since  the  output  values 
are  used  through  feedback  by  the  recursive  algorithm,  it  is 
possible  for  the  output  values  to  become  arbitrarily  large 
independent  of  the  input.   A  filter  of  this  kind  is  said  to 
be  unstable,  which  is  an  undesirable  condition.   Thus  we 
need  to  know  what  constraints  to  put  on  the  recursive  filter 
coefficients,  so  'that  the  filtering  operation  will  be  stable. 
This  problem  is  similar  to  the  stability  problem  for  the  one- 
dimensional  case,  except  that  the  added  dimension  inherently 
increases  the  complexity  of  the  analysis.   The  major  diffi- 
culty is  that  the  concept  of  zeros  and  poles  does  not  hold 
in  multi-dimensional  systems  so  that  simple  algebraic  systems 
extrapolation  of  one-dimensional  results  is  not  obvious. 
•  Stability  Theorems: 

Theorem  1  (Shank's  theorem):   A  causal  recursive  filter 
with  the  z-transform  H(ZwZ2)  =  A(z1 ,  z2)  /B  (z1 ,  z2)  ,  where  A 
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and  B  are  polynomials  in  z,  and  z2,  is  stable  if  and  only 
if  there  are  no  values  of  z,  and  z   such  that  B(z  ,z„)  =  0 
for  all  z.  _>  1  and  z«  >^  1.   (In  general  we  assume 

B(Z1'Z2}  =  b00  +  ^O2!"1  +  b01Z2"1  +  ^l2!*1^"1  ••" 
that  is  an  expression  in  increasing  powers  of  z.  and  z9.) 

In  other  words,  the  theorem  says  that  if  there  are  any 
values  (real  or  complex)  of  z,  and  z2  for  which  B(z,,z?) 
is  zero  for  z.  and  z2  simultaneously  greater  than  or  equal 
to  unity  in  magnitude,  then  the  filter  will  be  unstable. 

It  is  much  more  difficult  to  determine  the  stability  of 
two-dimensional  filters  than  of  one-dimensional  filters. 
In  the  one-dimensional  case,  it  is  only  necessary  to  locate 
finite  set  of  roots  in  the  z-plane.   In  order  to  mechanize 
the  theorem,  we  must,  in  general,  find  the  values  of  z.  and 
z_  for  which  B(z,,z~)  =  0.   There  is  no  technique  for  locating 
the  zeros  of  a  general  two-dimensional  polynomial  in  (z,,z_). 
Shanks,  Tretial  and  Justice,  in  their  paper  [11]  proposed  a 
technique  to  apply  theorem  1  to  determine  the  stability  of 
a  two-dimensional  recursive  filter,  H(z,,z2),  by  which  the 
unit  disk  in  the  z,  plane  must  be  mapped  into  the  z^  plane 
by  solving  the  implicit  two-variable  denominator  polynomial 
of  H(z,,z2).   A  necessary  and  sufficient  condition  for 
stability  is  determined  if  the  map  of  the  unit  disk 
d,  =    (z-,z_  _>  1)  does  not  overlap  the  unit  disk 
d2  E  (z?,z   >  1)  in  the  z   plane.   This  method  requires  an 
infinite  number  of  mappings  and,  therefore,  can  not  be 
applied  exactly. 
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Huang  [20]  simplified  Shank's  stability  theorem 
considerably  by  stating  the  stability  theorem  as  follows 
Theorem  2:   A  causal  filter  with  a  z-transform 


A(z,  ,z.J 
H(zwzn)   = 


'l'"2'      B(z1,z2) 


where  A  and  B  are  two-variable  polynomials  in  z,,z  ,  is 
stable  if  and  only  if: 

1)  the  map  of  the  unit  circle  6d   =  (z  , | z,  | =1 ) , 
according  to  B(z  ,z2)  =  0,  lies  outside  of  the  unit  disk 
d_  =  (z  ,|z  |_>1)  in  the  z  -plane,  and 

2)  no  point  in  the  unit  disk  d..  =  ( z ^  ,  |  z ,  | _>  1 ) 
maps  into  the  point  z   =  0  by  the  relation  B(z..,z2)  =  0. 

In  order  to  apply  this  theorem  as  a  stability  test,  we 
have  to  map  the  unit  disk  5d,  =  (z  , |z1 |=1)  into  the  z? 
plane  according  to  B(z,  ,z2)  =  0  or  z2  =  f  (z,)  |  i   i-,  and 
to  see  whether  the  contour  in  the  z^  plane  lies  inside  the 
unit  disk  d2  =  (z- ,  |  z2  |<_1)  .   Also,  it  is  necessary  to  solve 
B(z,,0)  =  0  to  see  whether  there  are  any  roots  with  magnitude 
greater  than  1. 

Although  theorem  2  is  much  simpler  than  using  Shank's 
original  theorem,  the  procedure  is  still  infinite  in  the 
sense  that,  in  principle,  we  have  to  map  the  unit  circle 
5d,  into  the  z? -plane. 

Ansell [21] ,  whose  main  contribution  is  to  couple  the  use 
of  a  Hermite  test  for  checking  stability  with  a  series  of 
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Sturm  tests  checking  positivity,  has  reduced  theorem  2. 
Although  Ans ell's  results  enable  us  to  test  stabilty  in  a 
finite  number  of  steps,  it  still  is,  unfortunately,  very 
tedious. 

Make  the  change  of  variables: 


and, 


1  -  Z2 

1  +  z2 


and  let 


E(p,  /P?) 
H(ZwZ0)   = 


l'~2 


fcp1,p2; 


where  E  and  F  are  polynomials  in  p,  and  p~.   We  can  restate 
Theorem  2  as  follows. 

Theorem  3  (Ansell's  theorem):   The  causal  recursive 
filter  H(z,,z2)  is  stable  if  and  only  if: 

1)  for  all  real  finite  w ,  the  complex  polynomial 
in  p„,  F(joj,p2)  has  no  zeros  in  ReP2  >  0,  and 

2)  the  real  polynomial  in  p-,,  F(p,,l)  has  no 
zeros  in  Rep,    0. 

This  theorem  is  essentially  the  same  as  Theorem  2,  but 
an  advantage  of  Theorem  3  is  that  condition  1)  can  be 
tested  using  standard  techniques  of  circuit  theory. 
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Anderson  and  Jury  [22]  modified  Huang's  method  by  out- 
lining a  procedure  which  replaces  the  bilinear  transforma- 
tions by  the  construction  of  a  Schur-Cohn  matrix  and  checking 
for  positivity  of  a  set  of  self-inversive  polynomials.   It 
also  replaces  the  Hermite  test  by  a  Schur-Cohn  matrix  test 
and  requires  a  series  of  Sturm  tests  or  equivalently ,  a 
system  of  tests  establishing  the  roots  distribution  of  a 
polynomial.   These  modifications  represent  a  substantial 
reduction  in  computations  as  compared  to  Huang's  method.   We 
will  discuss  stability  considerations  which  are  related  to 
our  design  technique  in  Chapter  IV  and  Chapter  V  in  detail. 

In  this  chapter,  basic  properties  of  two-dimensional 
digital  filters  have  been  reviewed  and  summarized  including 
semi-causality  to  establish  a  fundamental  theory  of  two- 
dimensional  digital  filters. 
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III.   MATHEMATICAL  BACKGROUND  FROM  ANALYTIC  PLANAR  GEOMETRY 

In  this  study,  we  propose  using  semi-infinite  planes  to 
approximate  a  given  frequency  response  of  a  two-dimensional 
analog  filters  by  extending  one-dimensional  semi-infinite 
straight  line  approximation  to  two-dimensional  cases.   For 
example,  as  we  will  see,  the  frequency  response  (dB  vs.  log- 
frequency)  of  the  separable  canonic  factor 

H(ja)1,j(u2)   =   cl+jUiTi)  (l+ja)2T2) 

can  be  approximated  by  a  collection  of  planes.   This  is 
discussed  in  Chapter  IV  in  detail. 

To  prepare  for  the  discussions  in  the  next  chapters, 
it  is  worthwhile  to  review  and  develop  useful  equations 
from  analytic  planar  geometry  which  are  going  to  be  used 
throughout  this  study  with  some  examples  of  how  to  use  them. 

The  properties  of  the  analytic  planar  geometry  can  be 
summarized  as  follows: 

1.   The  general  equations  for  a  plane  are  given  by 

Ax+By+Cz+D=0  (3.1) 

where  A,  B  and  C  are  the  direction  numbers  of  the  normal 
to  the  plane.   The  equation  of  a  plane  passing  through  the 
point  (xwy^z,)  and  perpendicular  to  the  line  which  has 
direction  numbers  of  A,  B  and  C  is  given  by: 
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A(x-x1)  +  B(y-y1)  +  CU-z^   =   0        (3.2) 

2.   The  direction  cosines  of  the  normal  to  the  plane 
are  given  by- 


cos  a   =   — ^ = T-TTo  (3.3a) 

(AZ  +BZ  +C  )1/Z 


cos  6   =   — ? ^ — TT7T  (3.3b) 

(Az  +B  +C  )±/Z 


cos  y   = 


{A2+B2+C2)l/2 


where  a,  8/  and  y  are  the  angles  that  the  normal  makes  with 
the  x,  y,  z  axes  respectively. 

The  relationship  between  direction  cosines  is  given  by 


2        2 

Cos  a  +  Cos  3  +  Cos  y   =   1  (3.3d) 


3.   For  a  line  joining  points  (x,  ,y^/Z-^)  and  (x2,y2>z2) 
the  equation  is  given  by 


Cos  a    _   Cos  3    _   Cos  y  ,^    ^, 

x2-xx    y2-y!     z2-z1 


4.   The  equation  of  a  straight  line  with  the  direction 
numbers  a,  b,  c  and  passing  through  the  point  (x^y^z^  is 
given  by 


x-Xl    y-y1    z-Zl 

a    ~     b     =     n —  (3.5) 


5.  The  angle  6,  between  two  intersecting  lines  with 
direction  numbers  a,  b,  c  and  a',  bf,  c'  respectively  is 
given  by: 

Cos  0   =   Cos  a  Cos  a1  +  Cos  6  Cos  3'  +  Cos  y   Cos  y1 

(3.6a) 
or 


Cos  9   =  aa'  +  bb'  +  cc'   -        (3.6b) 


Va2  +b2  +  c2  Va,2+b,2  +  c'2 


6.   For  parallel  lines,  the  relationship  between 
direction  numbers  is  given  by 


a'      b^     c*  U*/J 


7.  For  perpendicular  lines: 

aa'  +  bb'  +  cc '   =   0  (3.8) 

8.  The  direction  numbers  of  the  line  of  the  intersection 
between  two  given  intersecting  planes,  with  direction  numbers 
A,  B,  C  and  A' ,  B',  C  respectively,  is  given  by 
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a  = 


b   = 


B 


B*     C 


A 


C     A' 


(3.9a) 


(3.9b) 


c   = 


B 


A'     B' 


(3.9c) 


The  angle  6  between  these  two  planes  is  given  by 


Cos 


AA'  +  BB1  +  CC* 


V 


A2  +B' 


+  C   VA' 


(3.10) 


2     2     2 


9.   The  equation  for  a  plane  in  terms  of  its  slope  can 
be  calculated  as  follows: 

As  seen  from  Fig.  3.1,  the  slope,  m,  of  the  normal 
line  is  given  by 


m  =   tan 


h 

g 


^hFW 


(3.11) 


m  =   slope  of  the  plane 


M   =   tan (90°  -  9) 


=  -  2.     =   .  1 
h       m 


(3.12) 
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normal     to  olane 


normal      to     nlane 


plane 


Fig.  3.1.   A  General  Plane  Configuration 
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The  equation  of  the  plane  is  given  by  one  of  the 
following: 

A(x-x1)  +By  +  Cz   =   0  (3.13a) 

Ax  +  B(y-y,)  +  Cz   =   0  (3.13b) 


and 


Ax  +  By  +  C(z-z1)   =   0  (3.13c) 


A*l   =  By±      =      Cz1   =   -D  (3.13d) 


The  intersection  of  the  plane  with  the  z-x  plane 
is  given  by 


A(x  -x1)  +  Cz   =   0  (3.14) 


and  the  slope,  nu,  of  this  line  is 


»i   =  a§  =   -§  (3-14a) 


Similarly,  the  intersection  of  the  plane  with  the 
z-y  plane  is  given  by 

B(y  -y, )  +  Cz  =  0  (3.15) 
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and  the  slope  of  this  line,  nu ,    is 


m2   ~   "  dy~  "   "  C  (3.15a, 


Substituting  (3.15a)  and  (3.14a)  into  (3.13) 
yields  the  following  alternate  equations  for  the  plane. 

-m(x-x-j)  -  m2y  +  z  =  0  (3.16a) 

-m,x  -  m2(y-y2)  +  z   =   0  (3.16b) 

-itux  -  m2y  +  Cz  -  Zj)      =   0  (3.16c) 

The  direction  cosines  of  the  normal  in  terms  of 
m-,  /  m2  are 

-m, 

Cos  a   =   (3.17a) 


Vm^ 


2     2 


-m~ 
Cos  6   = (3.17b) 


v^7 


m22  +  l 


Cos  Y   = 


Vm-L2  +  m22  +1 


(3.17c) 


Substituting  (3.14a)  and  (3.15a)  into  (3.11)  and 
(3.12)  yields 


m  =  i-  (3-18a) 


m2 
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and 


M   =   -Vmi  +m2  (3.18b) 


2 


The  angle,  8,  between  two  intersecting  planes  from 
(3.10),  (3.14a)  and  (3.15a),  is  given  by 


m,m,  *  +  irunu  '  +  1  ,-  ,„, 

Cos  9   =   LJ: Li (3.19) 


ymj2  +  m22  +  1  Vmp+mp  +  1 


The  direction  numbers  of  the  line  of  intersection, 

from  (3.9),  (3.14a)  and  (3.15a)  are  given  by 

a   =   (m£  -  m2)  CC  (3.20a) 

b   =   (m£  -  m1)  CC*  (3.20b) 

c   =   (m|m2  -  m^ir^)  CC  *  (3.20c) 

The  direction  cosines  of  this  intersection  line 
are  given  by 


m2  ~  m2  (3.21a) 

Cos  a   =   


V2  2 

(ml  -  nu)  +  (m|  -  m,  )  +  (m-jnu  -  mlm^) 


2 


m'  -  m, 
Cos  6   = (3.21b) 


V2  2 

(mi  -nu)   +  ^mi  ~mi)  +  (m^m2  -m^ir^) 
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mlm2  "  m2ml ,,  ,.  , 

Cos  y   =   (3.21c) 


V^7 


2  2  2 

m2)  +  (m*  -m,)  +  (m|m  -  m'm, ) 


10.   When  two  planes  are  added  in  one  dimension  to  form 
a  new  plane,  the  following  results.   Given  two  planes 


-m,x  -  nuy  +  z  -  z,   =   0 


-m,'x  -  mly  +  z1  -  z'   =   0 


Their  sum  in  the  z-direction  with  z"  =  z*  +  z  is 
given  by 

-(m1+m')x  -  (m2+m2')y  +  z"  -  ^  +  z£)   =   0 

The  slopes  of  the  intercepts  with  the  zx  and  zy  planes 


m1"   =  -(m1  +  m|) 


m2"   =   -(m2  +  m2') 


The  slope  of  the  normal  is  given  by 

2  ,  ,_  ,  .,2,-1/2 


m 


"   =   [  (m,  +m[)   +  (m2  +m2)  ] 


(3.22 
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To  aid  in  comprehension,  two  examples  follow: 
Example  (1) :   Given  m,  (slope  in  z,x  plane),  m, 
(slope  in  x,y  plane)  ,  and  x,  (the  point  at  which  x-axis  goes 
through  the  plane).   Find  the  equation  of  the  plane. 
The  general  plane  equation  (3.1)  is  given  by 

Ax  +  By  +  Cz  +  D  =   0 

When  z  =  0 

B     D 
x       A  Y  "  A 


m3   =   If  =    "  I      (given) 


and 


When  y  =  0 


and 


x,   =  -  j  (given) 


A      D 
"   "  C  x  "  C 


■l  -  E  -  "I    (given) 


D 
Zl        C 
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If  we  substitute  these  into  the  general  equation,  the  result  is 


or 


A 

Ax  -  Am-,y  -  — z  -  Ax,   =   0 
3J    m,       1 


x  -  nuy  -  —  z  -  x,   =   0 
3    m,      1 

or 


z   =   nux  -  nunuy  -  ximi  (3.23) 

As  shown  here,  we  can  derive  the  plane  equation  easily 
from  the  given  specification.   After  finding  this  equation, 
we  can  find  the  canonic  factor  by  simply  letting 


x  =   log  go. 


y   =   log  a)- 

i  i  2 
z   =   10  log  |T| 

Alternately,  we  can  go  from  the  canonic  factor  to  equations. 
Thus  equation  (3.23)  can  be  written  as 

2 

10    log    |T(joo1,  joo2)  |         =      m1    log   oj1    -   m^    log    oj2    -    10    log   K 

(3.24) 
where 

10    log  K     =      ^ixi 
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This  result  assumes  that 

i  i2      wl 

|T(ju>,,ju>2)  |    = (3.24a; 

Taking  the  logarithm  of  (3.24a)  yields 

2 

10    log    |T|         =      10    p    log    u^    -    10    q    log    oj2    -    10    log   K 


(3.25) 


By  equating  (3.25)  and  (3.24),  we  get 


10  p  =  m1  (3.26a) 


10  q   =   m1m3  (3.26b] 


It  follows  from  (3.26a)  and  (3.26b)  that 


m3   =   q/p  . 


Thus  the  transfer  function  of  (3.24a)  has  a  logarithmic 
characteristic  as  sketched  in  Fig.  3.2. 

The  previous  result  can  be  extended  to  the  following 
transfer  function 


(1  +  jWl)p/2 

T(j<D,  /  ja>9)   =  —t-t? 7TT?  (3.27) 

1    Z       K±/Z(l  +  jw2)q/ 


By  substituting  (3.26a)  and  (3.26c)  into  (3.27),  then 
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Z=io  1^1T| 


■s  loovJ 


*=:  looi*/| 


Fig.  3.2.   The  logarithmic  characteristic  of 

P 


T(jgo1,  jo)2) 


jj- 


Kcu. 


20  log  |t|  =  m,  log  co,  -  m,  m   log  aj_  -  10  log  K 
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nu/20 

(1    +    jco,)     L 

T(>l'jw2)       ■      — Hu/20  (3*28) 

K1/Z(l   +    ja)2)     Z 


When   w,    <<    1   and   oj2    <<    1    ,    region   I 
T(joJl/jco2)       =      K~1/2 

and 

20    log    |T|       =      -    10    log   K 

When   co,    <<    1   and   oj-    >>    1*    region   II 


T(jUl,ja>2)       =      — ^-720 

K1/2(jaJ2)     A 


and 


20    log    |T|      =      -    10    log  K  -  nu    log   cj. 


When   oj2    <<    1   and   oj,    >>    1,    region   III 

111,720 

(j«2.) 

T(ju)1,jo)2)       =  -j/2 

K 


and 


20    log    |T|       =      -    10    log   K   +   m,    log 


0). 
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When  u>,  >>  1  and  oo-  >>  1,  region  IV 


nu/20 

j  OJ- 


T(ju),  ,jw~)   =   , 

1    2.  .  nu/20 


KX/2H   1  2' 


and 


20  log  |T|   =   -  10  logK  +  m1log  oj,  -  m,  log  to,. 

The  resulting  planar  approximation  for  the  logarithmic 
transfer  function  is  shown  in  Figure  3.3 
Example  2 ; 

In  this  example,  we  derive  planar  equations  from  the 
given  transfer  function. 

Given  a  separable  two-dimensional  transfer  function: 

T(jw,  ,  ju9)   = ■=         (3.29) 

(1  +  JU)1T1)P(1+  JU)2T2)H 

The  magnitude  square  of  the  transfer  function  (3.29)  is 
given  by 

K2 

|T(ja)1,ja)2)|  =      2 2~p" 2 2~q  <3-3°) 

1         A  (1  +  oJiZt1Z)P(1  +  ^2    T2    } 

Take  the  logarithm  of  both  sides  of  (3.30), 


10    log    |T|2      =    -10 p    log(l+co12x12)   -  10   q    log(l+u>2    x2    )   +  20    log   K 


(3.31) 


when   w,!,    <<    1   and   co2x2    <<    1,    the  equation    (3.31)    becomes 
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A  fcrtoUjITl 


t  7^ta>W 


IE 


Br 


3-  l°j^ 


>    TUlaj^l 


Fig.  3.3.   The  logarithmic  characteristic  of 

(1    +    jo))  Mi/20 

X         Z  K±/Z(l    +    jco2)il2/zu 

20    log     |t  I    =    -10  log   K   +  m,  log   u>,    -   m-    log   go- 
84 


2 

10  log  |T|    =   20  log  K   =   C 


where  C  is  some  constant. 

When   w^t^    >>    1   and   oj2t2    >>    1, 


or 


10  log  |T|    =   -20p  log  u^  -  20q  log  a,  T 


+  20  log  K 


i  i  2 
10  log  |T|    =   -20p  log  o^  -  20q  log  a>2  -  20p  log  t. 


Let 


-20q  log  x2  +  20  log 


K 


2 

z   =   10  log  |t| 


x   =   log  to. 


y  =   log  a). 


and  we  get  the  planar  equation, 


Z   =   - 


20px  -  20qy  -  20p  log  t,  -  20q  log  t2  +  20  log  K 


or 
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p  q 

T    T 

(-20p)x  +  (-20q)y  -  z   =   20  log  (-± ? ) 


K 


or 


n^x  +  m2y  -  z  -  v 


where 


m^   =   slope  of  the  line  (x,z)  plane 
=   -20  p  db/decade. 

m2   =   slope  of  the  line  (y,z)  plane 
=  -2  0q  db/decade. 

In  a  similar  manner  the  regions  (co,  t,  <<  1; 

CJ2T2  >>  "^  an<^  ^l1!  >>  ^'    (J°2T2  <<c  ^  can  ^e  invest:'-<?ate(i' 
The  planar  approximation  for  the  logarithmic  transfer  function 

is  sketched  in  Fig.  4.1  where  q  =  p  =  1. 

In  this  chapter  we  have  developed  some  useful  planar 

geometry  formulas  and  shown  how  they  can  be  used  to  determine 

a  planar  approximation  for  the  two-dimensional  logarithmic 

transfer  function  for  separable  cases.   Some  non-separable 

factors  are  discussed  in  the  next  chapter  in  which  experimental 

verification  of  this  approximation  technique  is  also  presented. 
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IV.   CANONIC  FACTORS  IN  CONTINUOUS  FREQUENCY  DOMAIN 

In  Chapter  III,  we  discussed  some  equations  of  planar  geo- 
metry and  showed  how  to  use  them  for  approximating  separable  two- 
dimensional  transfer  functions.   In  order  to  improve  the  scope 
of  the  proposed  design  technique,  we  now  consider  the  behavior 
of  other  separable  and  non-separable  factors  insofar  as 
their  planar  geometry  approximation  is  concerned.   This  will 
provide  insight  into  the  frequency  response  characteristics 
of  two-dimensional  systems  and  indicate  what  can  be  expected 
by  cascading  these  canonic  factors,  knowing  their  individual 
behavior.   The  logarithmic  characteristic  of  cascaded  factor 
is  simply  the  sum  of  their  individual  logarithmic  character- 
istics.  This  approach  parallels  the  semi-infinite  line 
approximation  (dB  vs.  log  w  )  technique  which  has  been  used 
so  successfully  in  one-dimensional  design. 

A.   SEPARABLE  FACTORS 

The  following  several  separable  canonic  factors  are  now 
considered. 

(1)  (1  +  jaj1T1)P(l  +  jt2u32)q  (4. a) 

i 

(2)  (1    +    jaJ1T1)P(l    +    ju)2T2)q(l    +    joo2T3)q  (4.b) 


(3)  [1   +    (jaJlT1)P]a      •  (4.c) 


where   p,    q   and   a   are    integers. 
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(1)   First,  as  an  example  of  the  technique,  consider 
the  following  analog  separable  transfer  function,  the  first 
separable  canonic  factor  with  p  =  q  =  -1; 

H(ju>  ,ju>0)   =   ,-,  7  ■ w-i  ■  • r  (4.1) 

1    2       (1  +  joj1t1)  (1  +  joo2t2) 

The  logarithm  of  magnitude-square  of  this  function 
is  given  by 

*>  *j  _  0    0 

10    log    |H(co1,oo2)  |         =      -10    logd+u^    t^2)    -    10    log(l+co2    t2    ) 

(4.2) 

Let  us  define  the  following  for  simplicity  of  notation. 

z   =   10  log  |H(oo1,o)2)  I    =   20  log  |H(o)1,a)2)  |  (dB) 

(4.3a) 
x  =   20  log  oo-l  (4.3b) 


y   =   20  log  w2 


(4.3c) 


Consider  the  regions  in  the  (gj,  ,u)2)  plane  (as  seen 
in  Figure  4.2)  as  follows: 

Region  I,  when  oo,x,  <  1  and  ^2x2  <  1, 

then  z  =  0   db  (4.4a) 

Region   II,    when   co^t  -^    >    1   and  co2t2    <    1 

z    =    -20    log    (w1t1)       db  (4.4b) 
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Region  III,  when  ^^l    <  1   and  W2T2  >  1 

z  =  -20  log(co2T2)  db  (4.4c) 


and 


Region  IV,  when  ou^  >  1  and  u   t~    >    1 


z  = 


=  -20  log  (oj1t1)  -  20  log  (oj2t2)  (4.4d) 

The  above  equations  (4.4a),  (4.4b),  (4.4c)  and  (4.4d) 
can  be  expressed  as 

Region  I,     z  =  0  (4.5a) 

Region  II,    z  =  -x  -  20  log  t-,  (4.5b) 

Region  III,   z  =  -y  -  20  log  t2  (4.5c) 

Region  IV,    z  =  -x  -  y  -  20  log  t2  -  20  log  x, 

(4.5d) 

These  equations  describe  four  intersecting  planes  in 
x,y,z  space  as  shown  in  Figure  4.1. 

Plane  I  is  the  horizontal  plane  at  zero  db.   Plane  II 
and  Plane  III  have  slopes  of  -20  db/decade  passing  through 
the  corner  lines  as  indicated  in  Figure  4.1.   These  lines 
are  analogous  to  the  break  point  frequencies  of  one-dimensional 
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Z  =   20  I  o  g  |H(  wr  w2  )  I 


og  w 


slope    —-20  d  B 


-   20dB 


Fig.    4.1.       Planar  Approximation    for 

1 


H  (co,  ,u>~) 


1+J0J1T1)    (1  +  J002T2) 


2       2 


20    log    I H  (oi1  ,co2)  |     =    -10    log    (1+w,     T, 

-10    log    (1+co22t22) 
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*2*1/^ 


Y  =  I  O 


g    w 


IT 


W 


I 


X   = 


°g  w. 


Fig.    4.2.      The   regions   of 


H(j<VJ0)2)      "      (1  +  J031T1)    (1  +  J0)2T2) 


2       2 


20    log    |H(o)1,o)2)  |       =    -10    log    (1+u^    Tj    ) 

2 
2    T2 


2       2 

-10    log    (1+oj„    t9    ) 


in    (log   u),  ,  log   co9)    plane 


91 


log-modulus  plots.   Several  properties  of  this  plane  are 
summarized.   Pertinent  equations  from  analytic  geometry  are 
given  in  Chapter  III. 

The  actual  frequency  response  of  this  separable 
canonic  factor  with  t,  =  t,  =  1  is  shown  in  Figure  4.3a 
(dB  vs.  log  w,,  log  oj2)  which  fits  the  planar  approximation , 
In  Fig.  4.3b/  the  contour  plot  of  the  frequency  response  is 
given,  which  verifies  the  regions  in  log  uj-,  /  log  u>-  plane 
shown  in  Fig.  4.2.   The  logarithmic  plots  versus  linear 
frequency  are  given  in  Fig.  4.3c  and  Fig.  4.3d.   It  is 
apparent  that  the  approximation  regions  are  well  defined 
when  logarithmic  frequency  is  used. 

The  frequency  response  of  this  separable  canonic 
factor  with  t,  =  1,  t-  =  2  is  given  in  Figures  4.4.   These 
results  again  confirm  the  logarithmic  approximation. 

(2)   Next  consider  the  canonic  form 

H(jMl,ja>2)   =   (1  +  j^tjj  (1  +  ju)2x2)  (1  + ja>2T3) 

where  we  have  two  separable  factors  of  oj2.   The  frequency 
response  of  this  factor  with  t,  *  1,    x2  =  2  and  x^  =  4  is 
shown  in  Figures  4.5.   As  seen  from  this  figure,  we  have 
-40  dB/decade  slope  in  the  direction  of  u2,  as  expected 
from  the  one-dimensional  case.   The  regions  are  shown  in 
the  contour  plot  with  dotted  line  (Fig.  4.5c) . 
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In  a  similar  manner  we  can  analyze  these  regions  as 
follows:   The  logarithm  of  magnitude  square  of  this  function 
is  given  by 


10  log  |H(u)1,u)2)  |2   =   -  10  log(l  +  oj12t12)  -  10  log  (1  +  w22x22) 


-  10  log(l +  w22t2) 


Consider  the  following  regions  in  (to,  , w2)  plane  (as 
seen  in  Fig.  4.5a)  as  follows: 

Region  I,  when  u-,t-,  <<  1  and  co2x2  <  ^2t3  <<  1 
then  z  =  0  dB. 

Region  II,  when  u,t,  >>  1  and  o)2t2  <  ^t-d  <k    1 

z   =   -  20  log  (oj-lt1)   dB 
Region  II,  when  co-,t-,  <<  1  and  oj2t2  <  ^2T3    <<;  1 

z   =   -20  log  (co2x2)  -  20  log(oj2x3)   dB 
Region  IV,  when  m,T,  >>  1,  ^2T3  >  C02T2  >:>  ^ 


z   =   -20  log  u,t,  -  20  log  ^2T2  "  20  l0g  W2T3 
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Region  V,    when   co^    <<    1,    ^^^    >>    1,    ^    T      >>o)2t2    >>    1 


Z      =      -20    logU2T3)    -    20    log(o,2T2) 


Region  VI,    when   ^^l    >>    ^'    W2T3    >    T2(J°2    >:>    1 


Z      =      -    20    log   a)1Tt   -    20    log    oo2t2    -    20    log   oj3t3 


where 


Z   =   20  log  |H(u)1,t2)  |  . 
(3)   Consider  the  separable  canonic  factor  now, 


H(a)1,o)2)   =   [1  +  (jw1T1)P]  a 


We  can  show  that  the  following  regions  occur  by 
using  the  notation  of  (4.3a)  and  (4.3b)  for  simplicity: 
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Fig.    4.3a.      The    frequency   response   of 

.  1 

H(jo)1,:o)2)       -       (i+jWi)  (l+jo)2) 

2 

z   =    2Q    log|H(o)1/o)2)  |    =    -10    logd+u^    ) 

(dB  vs.    log  bu,log   o>2) 


-10    log    (l+w2    ) 
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X-SCflLE=l .OOE+00  UNITS  INCH 
T-SCRLE=I . OOE+00  UNITS  INCH 


•ig.  4.3b.   The  contour  plot  of  H^,^)  =  (1+j   )  (i+j(0  ) 

z  =  20  log  |H|  =  -10  log  (l+o^2)  -  10  log  (l+io22) 
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Fig.  4.3c. 


The  frequency  response  of 
H(joJ1,jco2) 


(l+jai1)  (l+jo)2) 


z    =    20    log|H|    =    -10    log(l+oj12)    -10    logd+o^2) 
(dB   vs.    03,  ,0)-)     (linear    frequency) 
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10. 


>", 


a-, 


3. 


I 

T 


Fig.  4.3d.   The  contour  plot  of  H(oo..,go_)  =  . .  ,  . ,  .  =  ,  . r 

r  1   2  (1+]U)1)  (l+:co2) 

z  =  20  log|H|  =  -10  logd+u^2)  -  10  log(l+w22) 

(linear  frequency) 
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Fig.  4.4a.   The  frequency  response  of 
H(do)1,do)2) 


(l+joa1)  (l+j2co2) 


2  2 

=    20    log|H|    =    -10    log(l+oj1    )    -    10    log(l+4co2    ) 

(dB   vs.    log  co,  , log  co2) 
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X-SCRLE-1 .OOE+00  UNITS  INCH. 
T-SCflLE=l . OOE+00  UNITS  INCH. 


Fig.  4.4b. 


The  contour  plot  of  the  frequency  response  of 
H(o31w2)   = 


(l+ja)1)  (l+j2u)2) 


z    =    20    log|H|    =    -10    log(l+co12)    -    10    log(l+4to22) 


in    (log  to,  , log   to2)    plane 
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Fig.    4.4c.      The    frequency   response   of 

1 


H  ( j  co ,  ,  j  co  -, )       =  ,-,  ,  •  - 

J    1         2'  (l+jco,)  (1  +  j2co2) 

z   =    10    log|H|    =    -10  logd+co,2)    -    10    log(l+4co22 
(dB   vs.    co,  ,co_)        (linear    frequency   ccale) 
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Fig.  4.4d.   The  contour  plot  of  the  frequency  response  of 


H(jo)1,joj2)       = 


(1+w,     )  (l+4oo2    ) 
z   =    20    log|H|    =    -10    log(l+oo12)    -    10    log(l+4aj2    ) 
(dB   vs.    o)1,o)2)  (linear   frequency   scale) 
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Fig.  4.5a. 


The  region  of 
H(j(D1,ju)2)  = 


(l+jo)1T1)  (1+joo  t  )  (l+jco2x3) 
in  (log  a), /log  co„)  plane 
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Fig.  4.5b.   The  frequency  response  of 

1 
H(jo)1,jo)2)  -  (i+j0Jl)  (l+j2co2)  (l+34w2; 

z   =    20    log|H|    =    -10    log(l+0312)    -    10    log(l+4w2    ) 

2 
-10    log(l+16u>2    ) 

(dB   vs.    log   o)lflog  w2) 
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X-SCRLE=1 . OOE+00  UNITS  INCH. 
Y-SCRLE=I . OOE+00  UNITS  INCH. 


Fig.  4.5c.   The  contour  plot  of  the  frequency  response  of 
H(julfju2)  =  (i+jUl)  (l+j2oo2)  (l+34co2; 

2  2 

z  =  20  log|H|  =  -10  log(l+a)1  )  -  10  log  (1+4^  ) 

2 
-10  log(l+16oj2  ) 

(log  oo-L  vs.  log  co2) 
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Fig.  4.5d. 


z  = 


The  frequency  response  of 

H(joJ1,jOi2)  =  (1+j(Ui)  (i+j2co2)  (l+j4oo2) 

20    log|H|    =    -10    log(l+ca12)     -    10    log(l  +  2aj22) 

-10    log(l+16oo22) 
(dB   vs.    (ji-.,uiy)  (linear    frequency    scale) 
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10, 


0. - 


X- 


/«u 


Fig.  4.5e.   The  contour  plot  of  the  frequency  response  of 

.  1 

H(D(o1,30)2J    -     (1+jc0i)  (l+j2oa2)  (l  +  j4co2) 

?  2 

z   =    20    log|H|    =    -10    log    (1+w^)    -    10    log(l  +  4aj1    ) 

-10    log    (l+16oo22) 
(linear   frequency   scale) 
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Region  I,  when  go,  >  1/t,  , 

z   =   apx  +  20ap  log  t-,  dB 
Region  II,  when  to,  <  1/t-, 

z   =   0  dB 

These  separable  canonic  factors  are  identical  to 
the  canonic  factors  of  a  one-dimensional  transfer  function. 
They  are  just  an  extension  to  the  two-dimensional  case. 
Therefore,  they  need  not  be  discussed  further. 

B.   NON-SEPARABLE  FACTORS 

The  following  non-separable  canonic  factors  are  now 
considered. 

(1)  [1+  (jaj1x1)P(joj2T2)q]a  (4.d) 

(2)  [1  +  (ju1T1  +  joo2T2)P]a  (4.e) 

where   p,    q   and  a   are    integers. 

(1)       First,    consider   the   non-separable   canonic 
factor   of   the    form: 


H(oJ1,aj2)       =       [1  +  (jaj-LT1)P(joJ2T2)     ] 
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Region  I,  when 


P   q        1 

1      2  T  P  T  q 


z   =   0  dB 


Region  II,  when 


P   q 


1     2  T  p  ,  q 


Tl   T2 


z   =   apx  +  aqy  +  20ap  log  t,  +  20  ay  log  t-  dB 


The  cornerline  between  the  two  regions  is  given  by  the 
equation 


or 


p   q   p   q 

ul   w2   Tl   T2 


px  +  qy  +  20p  log  t,  +  20q  log  t2   =   0 


The  cornerline  is  sketched  in  Figure  4.6. 

The  slope  of  the  plane  in  region  II  is  given  by 


-42 


'p   +  q   .   Direction  numbers  of  the  normal  are  given 
by  -ap,  -aq  and  -1. 
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^  x  -  lo  g 
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y  -lo  g 


Fig.  4.6.   The  regions  of  non-separable  case  in 
log  to,  ,  log  co^  plane 
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Depending  upon  the  sign  of  p  and  q  and  the  magnitude  of 
t,  and  x2  several  different  orientations  of  the  corner line 
are  possible. 

Example  (1) 

Discuss  this  canonic  form  in  the  special  case  with 
p  =  q  =  1,  a=+l,  x,  =0.5,  and  x2  =  1.   The  transfer 
function  gets  the  following  form: 

H(o),,a)2)   =   [1  +  (.5jw,)  ( ja)2)  ] 

Take  the  logarithm  of  the  magnitude  square. 

i 2      |  i+2 

z      =      +20    log|l    -    .5(jj,co2| 
In   region   I,    w,t,    <<:    ^   an<^   aj?T2    <<:    "*" 


z      =      0      dB. 


In   region   II,    id,t,    >>    1   and   w2t2    >>    1 
z=x  +  y+20    log    (0.5) 


where 


x   =   20  log  co1,    y   =   20  log  cu2   and   z   =   20  log  |H(oj1,u2) 

The  actual  frequency  response  of  this  factor  is  shown  in 
Fig.  4.7  including  a  contour  plot. 

Ill 


Fig.  4.7a. 


The  frequency  response  of 

H(  joo1 ,  joo2)    =    1    +    ( ja>,)  ( jO.  5o32) 

z  =  20  log|H|  =  20  log  (1  -  0.5oj,oj2) 

(dB  vs.  log  go,,  log  oo^) 
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X-SCRLE=l;OOE+00  UNITS  INCH. 
T-SCRLE=1..0£iE  +  00  units    inch. 

Fig.  4.7b.   The  contour  plot  of  the  frequency  response  of 
H(jcoL,ja>2)  =  1  +  ( jco1)  (j0.5oj2) 
z  =  20  log|H|  =  20  log  (1  -  O.Su^u^) 
in  (log  co,  ,  log  ai2)  plane 
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When  a  =  -1  we  have  a  stability  problem  in  this  example 
The  denominator  of  this  transfer  function  becomes  zero 
when 


or 


1  —  0  .  5o)-j  0)2   =   0 


o)-|  cop   —   2 


This  is  the  singularity  of  the  transfer  function  in  (w,,^) 
plane-   It  can  be  called  a  singular  line  of  the  system.   It 
is  a  curve  in  the  (ui,,to2)  plane  to  be  compared  with  a  pole 
in  a  one-dimensional  transfer  function. 

Example  #2; 

In  this  example,  we  discuss  the  same  canonic  form 
with  different  coefficients,  i.e., 

p   =   l,       q   =   -1  ,       a   =   +l, 

x,   =   0.5,    x 2   =   1 

Then  the  transfer  gets  the  form  of 

*5ujl  +1 
H(Ul,«2)   =   (1  +  ~— > 
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Taking  the  logarithm  of  its  magnitude  squared 
|H(a),  ,u9)     =   ll  +  ±\Z 

.  5(jd, 

z   =   +20  log  ll  + 


w2 


.  5a).. 

In  Region  I  when  -  <<  1 

a)2 


z   =   0   dB. 

•  5ai, 

In  Region  II,  when >>  1 

w9 


z   =   20x  -  20y  +  20  log  (.5) 

where  x  =  log  u^,  y  =  log  u2  and  z  =  20  log  |h(u>,/u>2)  |.   The 
boundary  between  these  regions  is  given  by  the  equation 

.  5oi-|   =   oi~ 

We  also  have  a  stability  problem  for  this  canonic  factor  when 
u>2  s  0.  This  is  also  a  singularity  in  (u),  ,w2)  plane  in  the 
form  of  a  straight  line.   This  type  of  factor  is  unstable. 
(2)   Now,  consider  the  non-separable  canonic  factor 


H(o),,a>2)   =   [1  +  (ju),w,  +  ju^oO  ] 


We  have  two  regions  again,  as  follows: 
Region  I,  when  w,t,  +  w2t2  <  1 

z   =   0   dB 
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Region  II,  when  u),x,  +  w-x-  >  1 

z  =  ap  log  (oa1x1  +  aj2T2)  • 

The  separation  line  between  the  two  regions  is 
given  by  the  equation: 


U1T1  +  0)2T2   =   1 


or 


log  (10,  x,  +  o)2t2)   =   0  .  (4.6) 

The  locus  of  points  of  (4.6)  is  shown  in  Figure  4.8, 
For  region  II,  we  have  to  think  of  two  conditions, 
as  follows: 

a)  When  w,x,  >  oo2x2 

z   =   apx  +  ap20  log  x, 

b)  When  oj,t,  <  ^2x2 

z   =   aqy  +  ap20  log  x2  . 

Thus,  region  II  is  split  into  two  areas  with  the  separation 
line  given  by 


W1T1   =   W2T2  ' 
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Fig.  4.8.   Locus  points  of  logC^T^  +  w2x2)  -  0 
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Frequency,  response  of  this  non-separable  transfer  function 
is  shown  in  Fig.  4.9  with  p  -  1,  a  =  -1,  x  =  x   =  l. 

C.   CONCLUSIONS 

In  this  chapter,  we  discussed  the  properties  of  the 
separable  and  non-separable  two-dimensional  properties  of 
canonic  factors  in  the  analog  domain  (co,,^  plane).   We  can 
conclude  that  these  canonic  factors  can  be  approximated  by 
planes  in  the  log  magnitude  versus  log  co,  and  log  oj^  space. 
As  a  design  technique  we  can  obtain  any  given  frequency 
specification  by  proper  combinations  of  planar  approximation 
of  canonic  factors.   It  is  obvious  that  only  certain  special 
cases  and  symmetries  can  be  obtained  using  the  factors  dis- 
cussed in  this  thesis.   For  example,  a  wedge-shaped  character- 
istic in  frequency  domain  could  not  be  synthesized  by  means 
of  the  factors  we  have  considered. 

The  other  result  of  this  study  is  the  stability  problem 
which  arises  in  certain  transfer  functions,  that  is, the 
singularity  curve  or  line.   For  a  separable  case,  we  don't 
have  a  stability  problem;  but  in  non-separable  cases  we  may 
have  stability  problems. 
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Fig.  4.9a. 


The  frequency  response  of 
J  1  J  2     1  +  joo,  +  jog. 


z  =  20  log|H|  =  10  log  1 1+ (oj, +oj2) 
(dB  vs.  log  to,  /log  co2) 


-,  2 


119 


X-SCRLE=1  . OOE  +  00  UNITS  INCH 
T-SCflLE-1  . OOE+00  UNITS  INCH 


Fig.  4.9b. 


The  contour  plot  of  the  frequency  response  of 
H(ju>,  ,jw 0)    =      . —-. 

z  =  20  log|H|  =  10  log  |l  +  (u^+u^)2! 
in  (log  co,  ,  log  co-)  plane 
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Fig.  4.9c.   The  frequency  response  of 

1 
H(3oa1,:a)2)  = 


1  +  ju). 


JO). 


z  =  20  log|H|  =  10  log  |l  +  (u^+a^)  | 
(dB  vs.  a)1,oo2)    (linear  frequency  scale) 
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v-   TWO-DIMENSIONAL  DIGITAL  FILTER  DESIGN  TECHNIQUE 

A.   INTRODUCTION 

The  design  of  two  and  multi-dimensional  recursive  digital 
filters  is  difficult  due  to  the  absence  of  the  Fundamental 
Theorem  of  Algebra  in  two  or  more  dimensions  [15].   That  is, 
polynomials  in  two  variables  may  not,  in  general,  be  factored 
into  lower  order  polynomials.   For  this  reason,  several  design 
procedures  for  two-dimensional  recursive  filters  have  been 
proposed  in  the  literature.   We  then  consider  the  use  of  the 
double  bilinear  z-transfer  with  the  semi-infinite  plane 
approximation  in  the  previous  chapters.   These  are  briefly 
discussed  in  the  following. 

Shank  et  al  suggested  three  design  techniques  in  their 
paper  [11] . 

1.  The  use  of  a  series  of  stability-preserving  transforma- 
tions to  produce  a  stable  two-dimensional  recursive  discrete 
filters  from  a  stable  one-dimensional  continuous  filter. 

2.  A  space  domain  design  technique  in  which  the  unit 
sample  response  of  the  filter  being  designed  is  made  to 
approximate  a  given  desired  unit  sample  response. 

3.  A  stabilization  technique  in  which  the  denominator 
polynomial  is  replaced  by  its  double  planar  least-squares 
inverse.   Planar  least-squares  can  be  defined  as  follows: 

We  have  some  given  array  C.   We  would  like  to  find 
an  array  P  such  that  C  convolved  with  P  approximates  the 
unit  pulse  array  U.   That  is, 

C  *  P  «  U 
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where  the  symbol  *  denotes  two-dimensional  convolution.   In 
general  it  will  not  be  possible  to  make  C  *  P  exactly  equal 
to  U.   In  actuality,  C  *  P  will  be  equal  to  some  other 
array,  G.   If  we  choose  P  such  that  the  sum  of  the  squares 
of  the  elements  of  U-G  is  minimized,  P  is  called  a  planar 
least  square  inverse  of  C.   Planar  least  square  inverse  of 
P  is  called  double  planar  least  square  inverse  [11] .   This 

new  polynomial  is  conjuctured  to  have  no  zeros  on 

-2  A  ,    ,     ,   , 

u   -   {z^z^:  |  ZjJ  ^_,  |z2|  >_  1}  and  to  have  a  frequency 

response  which  closely  approximates  that  of  the  original 
transfer  function. 

The  first  technique  is  known  as  a  filter  design  with 
rotating  axes  and  these  filters  are  called  rotated  filters. 
This  design  procedure  consists  in  considering  a  transfer 
function  in  the  analog  frequency  domain,  rotating  the  axes, 
and  applying  the  bilinear  transformation  to  obtain  a  two- 
dimensional  recursive  filter.   In  general,  the  starting 
point  is  a  separable  filter  and  stability  of  the  resulting 
filter  is  not  guaranteed  even  if  the  prototype  filter  is 
stable.   This  technique  was  used  by  Costa  and  Venetsanopoulus 
[23]  to  design  low-pass  filters  having  circular  symmetry, 
and  also  they  proved  that  rotated  filters  are  stable  if  the 
angle  of  rotation  is  between  270°  and  360°.   A  technique 
similar  in  character  but  rotated  axes  in  digital  frequency 
domain  ((z,,z2)  domain)  has  been  given  by  Chang  and  Aggarwal 
[24].   Hirona  and  Aggarwal  [25]  have  accomplished  the  rota- 
tion through  the  introduction  of  rational  powers  of  z^  and 
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z2,  and  the  interpolation  of  input  and  output  sequences. 
They  have  proved  that,  in  this  technique,  the  resulting 
filter  is  stable  if  the  prototype  filter  is  stable.   Further, 
the  distortion,  a  consequence  of  the  bilinear  transformation, 
is  absent  in  the  resulting  filter  when  compared  with  Shank's 
et  al. 's  filter. 

The  authors  of  [11]  were  unable  to  prove  that  the  new 
denominator  polynomial  produced  by  the  third  technique, 
which  is  an  extension  of  the  one-dimensional  case,  was 
stable  in  the  two-dimensional  case.   Nonetheless,  several 
hundred  filters  have  been  designed  using  this  technique 
without  a  counter  example  appearing  [18];  furthermore,  Jury, 
Kolavennu  and  Anderson  [26]  have  shown  it  to  be  valid  for 
certain  low  order  cases.   However,  Genin  and  Kamp  [27]  have 
recently  shown  that  counter-examples  do,  in  fact,  exist. 

Spectral  factorization  is  another  technique  which  can  be 
used  to  produce  a  stable  filter  from  an  unstable  one.   By 
spectral  factorization  we  mean  a  complex  map  that  carries  a 
stable  rational  transfer  function  into  another  stable  transfer 
function  exhibiting  a  different  frequency  response,  at  the 
same  time  maintaining  some  desirable  characteristic  [29] . 
Pendergrass,  Mitra  and  Jury  [28]  used  stability-preserving 
spectral  transformations  to  modify  the  frequency  response  of 
quarter-plane  filters. 

Mitra  and  Chakrabarti  [29]  extended  this  work.   Read 
and  Treitel  [30]  suggested  a  two-dimensional  Hilbert  trans- 
form to  perform  the  spectral  factorization,  while  Pister 
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[31]  and  Ekstrom  and  Woods  [32]  made  use  of  the  two- 
dimensional  cepstrum.   Ekstrom  and  Woods  also  suggested  ways 
that  the  factors  could  be  truncated  and  smoothed  to  preserve 
stability.   They  also  showed  that  the  usual  quarter-plane 
filter  is  not  the  most  general  two-dimensional  recursive 
filter.   It  is  instead  possible  to  have  a  recursive  filter 
which  in  the  limit  has  a  unit  sample  response  with  support 
on  an  asymmetric  half-plane,  and  also  they  developed 
stability  conditions  for  this  new  class  of  filters. 

B.   TWO-DIMENSIONAL  DIGITAL  FILTER  DESIGN  PROCEDURE  BY 
DOUBLE  Z-TRANSFORM 

1.   Double  Z-transform 

In  the  one-dimensional  recursive  digital  filter 
design,  the  bilinear  z-transform  is  often  applied  to  an 
analog  filter  transfer  function.   We  proposed  that  this  one- 
dimensional  design  approach  can  be  extended  to  the  two- 
dimensional  case.   Before  going  into  the  design  procedure, 
two-dimensional  recursive  filter  design  via  the  double 
bilinear  z-transform,  it  is  better  to  prove  that  we  can 
use  double  z-transform  for  two-dimensional  systems. 

Each  of  the  canonic  factors,  proposed  in  Chapter 
IV,  either  separable  or  non-separable,  can  be  converted 
into  discrete  equations  by  the  double  z-transform,  as 
follows: 

2  1  -  ZI' 
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2  X  "  z2 
1  +  z2 


The  single  z-transform  is  equivalent  to  trapezoidal 
integration  which  can  be  shown  as  follows: 

X(s)   =   S  Y(s)  (5.2) 


.<*>    =  &J§i. 


or 

n 

Yn  ~  yn_x   =     /  x(t)  dt 
n-1 


?(x   +  x   .)  (5.4) 

2   n    n-1 


where  T  is  the  sampling  period. 

We  write  (5.4)  in  the  z  transform  domain  as 


(1  -  z  1)Y(z)   =   |d  +  z  X)X(z) 


Hence 


X(z)   =   §(1  "  \[)     Y(z)  (5.5) 

1  +  z 
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Comparing  (5.5)  with  (5.2),  it  is  seen  that 

2  1  -  z"1 
S   =   J(i ^-T)  (5.6) 

1  +  2 

Equation  (5.6)  is  known  as  a  trapezoidal  bilinear  z -trans form, 

We  can  drive  the  double  z-transform  in  a  similar 
way,  as  follows: 


or 


Then 


,2 
W(x,y)   =  33^  Z(x,y)  (5.7a) 


W(s,,s2)   =   sis2  z(s]_'s2^  (5.7b) 


/  /d   Z(x,y)   =   /  /W(.x,y)  dx,  dy       (5.8) 


by  a  discrete  approximation,  the  left-hand  side  of  (5.8) 
becomes  (see  Fig.  5.1) 


/  /d2  Z(x,y)   =   [Z(xn,yn  "  z(xn'yn-l)] 

"  ^V-l'V  "  Z(xn-l'Yn-l)] 
-1   -1  .  -1  -1, 


(1  -  zx  -  z2  +  z1  z2  )Z(.z1,z2) 


(5.9) 
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and  the  right-hand  side  of  (5.8)  is  approximated  as  follows: 


//W(x,y)dxdy  =  A^£[w(Xn,yn)  +  W(xn,yn-1) 


+  "(Vl'V  +  W(xn-l'*n-l)] 


AXAVr,  -1    ,       -1    ,       "I   -1,TT/  V 

=  — -^- [1  +  z,   +  z2  +    zi    z2  ]W(z1,z2) 


(5.10) 


By  equation  (5.10)  and  (5.9)  we  have 


W(z1,z2)      1  -  z"1  -  z"1  +  z~  z2 


Z(z,,z0)      .     -1     -1     -1  -1    AxAy 
1  2  1  +  z,   +  z2   +  z,  z2       2 


(5.11) 


or 


-1     i     -1 

2  1  "  Zl    2  X  ~  Z2 

=   ^(T—*^)  ^(TT7^)  (5*12) 

1  +  z,     J  1  +  z2 


By  comparing  (5.1a)  and  (5.2b)  it  is  seen  that 

2  1   "  zl 

1  "  Z2 

Q.E.D. 

Thus,  this  is  called  a  double  z-transform.   The 
double  z-transform  seems  to  be  a  reasonable  substitution 
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y(t) 


x  (O 


—       n-1  n 


►    t 


(a)      One-dimensional   case 


Z(*,y) 


w(x,y) 


>W(x,y) 


(b)   Two-dimensional  case 


Fig.  5-1-   Two-dimensional  extension  of  the 
Bilinear  Transformation 
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for  each  canonic  factor.   However,  each  factor  must  be 
checked  for  stability,  which  is  discussed  in  detail  later  in 
this  chapter.   If  each  transformed  factor  is  BIBO  stable  the 
overall  transfer  function  will  be  BIBO  stable.   But  sometimes  the 
double  z-transform  may  in  certain  cases  lead  to  unstable 
solutions  [33] .   Finally,  it  is  noted  that  the  double  z-transform 
results  in  a  quarter-plane  symmetry  which  is  a  severe  limitation. 
2 .   Design  Procedure 

The  two-dimensional  recursive  digital  filter  design 
investigated  in  this  thesis  is  the  double  bilinear  z-transform. 
Proceed  as  follows: 

1.   Determine  a  two-dimensional  analog  transfer 
function, 

A(s1  ,s9) 

hi«i'-2>    =  niTTifr  (5-12> 

where  A(s,,s2)  and  B(s,,s2)  are  mutually  prime  two-dimensional 
polynomials  and  B(s1,s2)  satisfies 

B(s1,s2)  f      0 

(5.13) 

V(s1,s2)  z    r   =   {(Sl,s2):   Re  (s1)  >_0,Re(s2)  >_  0 } 

in  such  a  way  that,  the  frequency  response  of  G(s1,s2) 
meets  the  given  specifications  (for  example  as  discussed 
in  Chapters  III  and  IV) . 
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2.   Apply  the  double  bilinear  z-transform 

2  1   "  zl1 

51  =   A{ ^T}  (5.14a) 

1  +  z, 

2  1   "  zl 

52  =   -(- ij)  (5.14b) 

1  +  z2 

where  A  is  the  sampling  period,  which  is  assumed  the  same  in  both 
directions.   This  yields  the  two-dimensional  digital  transfer 
function 


2  1  -z,      1-2,  C(z,,z  ) 

H(f( TI)'f( =T)}   =  T^VZ2]    ~      D(z  z    )  (5'15) 

A  1  +  z2X   A  1  +  Z21  X   2     D(Z;l,z2j 


where  C(z,,z2)  and  D(z,,z2)  are  mutually  prime  two-dimensional 
polynomials  in  z,  and  z2.   (5.15)  will  have  the  desired 
frequency  response  with  proper  sampling  time  without 
aliasing  effect.   Since 


,  1  -  z?1 

si  =  !< *r 

1     A  1  +  z. 


maps  the  closed  unit  disc  onto  the  closed  right  half 
plane,  one  would  expect  (5.15)  to  satisfy  the  following 
sufficient  stability  conditions  [11],  [20]. 


D(z1,z2)  ?      0 
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\AW  £  u2   -  {(.zvz2):    |Zl|  >_  1, 


z2|  >_  1}     (5.16) 


Unfortunately,  however,  the  conditions  (5.13)  and  (5.16) 
as  we  will  show  later  in  this  chapter,  are  not  always 
satisfied. 

In  summary,  we  have  two  conditions  to  be  satisfied 
for  stability. 


I.   B(slfs2)  ?      o 


\/(s1,s2)  E  r   =   {(s1,s2):  Res1  >_  0,Res2  >  0}  (5.17) 


II.   D(21,z2)  ?      0 


Y(z1,z2)  £  u   =   {(z1,z2):  Iz-J  >  l/|z2|  >_  1} 

(5.18) 
where  V  stands  for  "for  each",  and  u   is  the  closed  unit  bidisc 

C.   DIGITAL  CANONIC  FACTORS  AND  STABILITY  CONSIDERATIONS 
1 .   Introduction 

Certain  properties  of  two  variable  polynomials  and 
rational  functions  which  will  be  needed  later  will  now 
be  discussed.   It  is  well  known  that  a  two  variable  polynomial 
is  not,  in  general,  factorable  into  first-order  polynomials; 
but,  a  two  variable  polynomial  can  be  factored  into  irreduci- 
ble factors  which  are  themselves  two  variable  polynomials 
but  which  cannot  be  further  factored.   Of  course  a  given 
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polynomial  may  itself  be  irreducible.   These  irreducible 
polynomials  are  unique  up  to  multiplicative  constants  [33] . 
Two  polynomials  which  have  no  irreducible  factors  in  common 
are  said  to  be  mutually  prime.   Consider  two  variable 
rational  functions,  in  general, 


P (z, ,z9 
G(z, ,z0)   = 


l'"2'      Q(z1,z2) 


where  P(z^,z2)  and  Q(z,,z2)  are  mutually  prime.   Using 
terminology  of  [34],  a  point  (zj,zl)  making  Q(z|,z2)  =  0 
but  P(z',z')  ^  0  is  called  a  pole  or  a  nonessential  singu- 
larity of  the  first  kind  (such  a  point  is  analogous  to  a 
pole  in  the  one  variable  case).   A  point  (z£,z£)  making 
P(zV,z~)  =  Q(zV,z2)  =  0  is  called  a  nonessential  singu- 
larity of  the  second  kind  (such  a  point  has  no  analogue  in 
the  one  variable  case).   Clearly,  if  (zi,zl)  is  a  pole, 
G(z,',z')  =  °°,  but  if  (zV,z")  is  a  nonessential  singularity 
of  the  second  kind,  then  G(zV,z")  is  undefined. 
2.   Separable  Factors 

The  most  common  canonic  separable  factor  in  Chapter 
IV  is 


H(<jj,,u)0)   —   "7t — —j, \  1 1  .  a r        (5.19a) 

1   2       (1  +  Joj,  x,)  (1  +  3oj2t2) 


assuming  x   =  x_  =  1.   Substituting 


jw-L   -   S1  , 
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ju>2      =      S2 


We   get 


H(si's2»     "     (TTT[flTTip  <5-19b> 

Substituting  equations  (5.1a  and  5.1b)  into  (5.19b), 
performs  the  double  z-transformation.   Thus, 

8(1  +  z"1)  (1  +  z"1) 

T(zrz2)   =   i- £_  (5.20) 

(a  +  z1    )  (a  +  z2  ) 


where 


A2 


(A  -2)2 


a   = 


a  +  2 
a  -  2 


A   =   sampling  period 

This  separable  canonic  form  obviously  satisfies  the  condition 
of  (5.17).   We  should  be  careful  about  the  second  condition 
of  (5.18) . 

When  z,  =  z2  =  -1,  we  have  a  nonessential  singularity 
of  the  second  kind  if  a  =  1,  but  a  =  1  never  occurs. 
Therefore,  there  is  no  nonessential  singularity  of  the  second 
kind.   Thus,  this  separable  case  is  stable  by  satisfying 
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two  essential  necessities  (5.17)  and  (5.18).   The  frequency 
response  of  this  separable  case  is  plotted  in  Fig.  5.2,  with 
sampling  time  A  =  0.08  which  fits  that  of  the  analog  frequency 
response. 

The  double  z-transform  for  equation  (5.19a)  is 


A2  (l  +  z11)  (l  +  z^1) 


T(Z1/Z2)   =   Ta-2T;l)  (A-2t2)    A+2T;l     ~   aT2~^     ~ 

(A^2T7+Z1  )(A^2TJ+Z2  > 

(5.21) 

and  this  general  case  for  the  two-dimensional  separable 

factor  is  stable  by  satisfying  the  stability  necessities 

of  (5.17)  and  (5.18) 

The  second  analog  canonic  factor  in  Chapter  IV  is 


H(jUlrJ0.2)  (1  +  jo)1T1)  (1  +  JW2T2)  (1  +  j  w2T3)     ' 

When  we  perform  the  double  bilinear  transform  for 
this,  we  get 


3      -1       -12 
AJ(1  +  z/)  (1  +  z9xr 

T(zn,z„)  = 


1'  2'  A+2t,     ,   A+2T-     ,   A+2t-.     , 

(A-2x1)(A-2x2)(A-2x3)[^+z-i][I=2^+z2  HI=ni+«;  ] 

(5.22) 

The  frequency  response  of  (5.22)  for  t,  =  1,  t2  =  2,  t^  =  4 
is  shown  in  Fig.  5.3  which  almost  fits  the  analog  case. 
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Fig.  5.2a.   The  frequency  response  of 


T(z1,z2)  = 


where  3  = 


e(l+z1"1) (l+z2  1) 

(a+z^  )  (a+z~  ) 

A2  A+2 
~  ,  a  = 


(A-2)" 

(dB  vs.  log  a),  ,  log  u}~) 


A-2 


=  0.08 
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Fig.  5.2b.   The  frequency  response  of 


T(z1,z2)  = 


6(l+z1~1) (l+z2  1) 
(a+z1"1) (a+z2-1) 


3  = 


A' 


(A-2) 


2  ' 


a  =  |i|  ,  A  =  0.08 


(dB  vs.  a).  ,0)-)    (linear  frequency  scale) 
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IT). 


0.  -, 


*4 


10.-, 


Fig.  5.2c. 


The  contour  plot  of  the  frequency  response  of 
3(l+z1"1)  (l+z2_1) 


T(zx,z2) 


6  = 


2  ' 


(a+z1~  ) (a+z2   ) 
a  =  f±§,  A=  O.OB 


where 


(A-2) 
(linear  frequency  scale) 
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We  can  summarize  the  m  xn  analog  separable  factor  as 


M              N 
H(ju>  ,ju,  )   =   [  IT   ,    .-*• ]  [  n   ,  ,  / ]        (5.23) 


After  performing   double   z-transform,     (5.23)    becomes 


M  -1  N  -1 

M  a    (1  +  z,    )  N  A    (1  +  z,) 

T(z,,z    )       =       [    n      — =i ][    n      — | ] 

(A-2t.)  (-— ii+zT1)       3    L     (A-2x.)  (_-^J-+z:1) 
1       A-2x  -1  3        A-2t  ■  2. 

15.24) 

Equations  (5.23)  and  (5.24)  satisfy  the  conditions 
of  (5.17)  and  (5.18).   As  a  result,  a  separable  filter 
results  in  separable  filters  in  digital  domain  and  they  are 
stable. 

3.   Non-separable  Factors 

The  first  non-separable  factor  in  Chapter  IV  is 
4.d  with  a  =  -1,  p  =  q  =  t,  =  t2  =  1  as  given  below 


1 


H  (  joj,  t  j  oj0)   -  j— -  ■ 

-1     1    2        1  +  3oj-i](jo~ 


(5.25a) 


or 


«<VS2>  ■  TTT^  (5-25b: 
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Fig.  5.3a.   The  frequency  response 

A3(l+z1"1)  d+z^1) 
T(z1/z2)  = 


(A-2)  (A-4)  (A-S,  (f±|  .z,-1)  (f±|  .zf1)  <£±f  +z2-h 


where  A  =  0.0  8 
(dB  vs.  log  co,  ,  log  u^) 
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Fig.  5.3b.   The  frequency  response  of 


T(2l,z2)  = 


A3(l+z1"1)  U+z^1) 


(A-2) (A-4) (A-8) (|±|  +z1"1) (|±i  +zx  X) (f±f  +  z2  1) 


where  A  =  0.08 


(dB  vs.  u.,ii)J 
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!S.u 


Fig.  5.3c.   The  contour  plot  of  the  frequency  response  of 


T(z1,z2)  = 


A3(l  +  z1"1)  (l+z^1) 


(A-2)  (A-4)  (A-8)  (|±|  +z1-1)  (|±J  +a1"1)  <£±§  -hz^1) 


where  A  =  0.0  8 
(linear  frequency  scale) 
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and  by  the  double  2 -trans form,  we  get 


A2(l  +Z,1)  (1  +  Z'1) 

T(Z1'Z2)   =  — 2 =T TT T Tt-ZT   (5'26) 

(A  +4)  (A  -4)  (z  L  +  z  l)   +  (A  +4)z  1z  X 


In  the  analog  domain  (5.25b)  does  not  satisfy  the 
condition  (5.17).   We  have  a  singularity  at 


Ulu)2  =   1 


For  the  second  condition,  we  have  a  nonessential 
singularity  of  the  second  kind  at  z,  =  -1,  z_  =  -1.   This 
does  not  tell  us  whether  it  is  stable  or  not,   When  we 
check  the  frequency  response,  it  is  seen  that  this  factor 
is  unstable. 

The  second  non-separable  canonic  factor  can  be 
derived  from  (4.d)  with  a  =  q  =  -1  and  T-,=T2  =  las 
follows. 


H(ju)i  ,jui0)   =   3 (5.27a) 

1  + 


2  tu. 


or 


H<si<s2>    -   — 4r  =   sr^r  <5-27b> 

1  +  IT 
s2 
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Then,  (5.27b)  becomes  after  performing  the  double  z-transform, 
as  follows: 

(l  +  zT1)  (1  -z'1) 

T(z  ,z  )   = _     2  (5.28) 

2U-«1r«2I) 


This  canonic  factor  does  not  satisfy  either 
stability  condition.   Hence  it  is  unstable. 

The  other  series  of  canonic  factors  can  be  derived 
from  (4.e) . 

For  example,  with  a  =   -1,  t1  =  t-  «  p  »  1,  we 

get 


HtjoiwjwJ   =   -,  ,  -1  ,  ■ (5.29a; 

-L    ^       1  +  3  co,  +3  u^ 


or 


H(S1'S2»   =   l  +  s^  +  s2  (5-29b> 


or  in  the  digital  domain 


A(l  +  z"1)  (1  +  z"1) 


T(z,,z  )   =   — —, _        (5.30) 

(A+4)  +A(zxx  +  z2  )  +  (A-4)z1  z2 


This  canonic  factor  satisfies  the  first  condition  but 
it  has  a  nonessential  singularity  of  the  second  kind,  at 
z,  =  z~  =  -1  and  so  we  should  check  the  frequency 
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response.   We  see  that  there  is  a  singularity  at 

(u)1,o)2)  =  (it,  it).   Thus  this  non-separable  canonic  factor 

is  unstable. 

4 .   Conclusions 

All  of  the  canonic  factors  discussed  in  this 
chapter  are  collected  in  Table  I.   By  comparing  the  results 
shown  in  this  table,  we  can  say  that  separable  stable  canonic 
factors  result  in  stable  separable  filters.   As  seen  from 
Table  I,  the  separable  transfer  functions  1,  2  and  3  which 
start  with  stable  separable  ones  in  the  two-dimensional 
continuous  frequency  domain  end  up  as  stable  separable  trans- 
fer functions  in  the  two-dimensional  discrete  domain.   This 
is  also  true  for  the  mxn  separable  canonic  factor  4. 

Non-separable  canonic  factors  have  very  different 
characteristics.   The  unstable  non-separable  case  ends  up 
with  a  non-separable  case  in  the  discrete  domain,  which  is 
also  unstable.   Although  we  could  not  prove  this  result 
generally,  we  also  could  not  find  any  counter  example  to 
exist.   The  transfer  functions  of  5,  and  6  start  in  the 
unstable  case,  and  end  up  with  the  unstable  case  in  the 
discrete  domain.   The  transfer  functions  of  7  starts  with 
the  stable  case  in  the  continuous  frequency  domain,  ends  up 
with  the  unstable  case  in  the  discrete  domain.   Thus  we 
have  to  be  very  careful  with  the  double  z-transform  since  it 
does  not  yield  stable  results  every  time.   It  may  yield 
unstable  results  even  if  it  starts  from  the  stable  case. 
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VI.   CONCLUSIONS 

A  technique  for  the  use  of  semi-infinite  planes  to 
approximate  the  log  magnitude  versus  log  frequency  charac- 
teristics of  two-dimensional  filters  is  presented.   This 
technique  is  applied  to  quarter  plane  filters  and  works 
well  for  separable  transfer  functions.   Other  non- 
separable  canonic  (basic)  transfer  functions  are  also 
studied  in  terms  of  their  planar  approximation. 

The  following  separable  and  nonseparable  canonic  factors 
have  been  studied. 

(1)  (1  +  jo)1x1)p(l  +  jw2T2)q 

(2)  (1  +  oo1x1)P(l  +  jco2T2)q(l  +  joj2T3)q' 

(3)  [1  +  (jco1t1)P]a 


q,a 


(4)    [1  +  (joo1T1)^(l  +  ja>2T2)*] 


p,a 


(5)   [1  +  (ja)1T1  +  jco2T2)^] 


The  properties  and  stability  of  these  canonic  factors 
are  investigated  first  in  the  continuous  variable  domain 
and  then  after  the  double  bilinear  z-transform  is  applied 
to  extend  them  to  two-dimensional  recursive  digital  filters 
These  are  all  summarized  in  Table  I  in  Chapter  V.   An 
important  result  is  that  although  the  double  bilinear 
z-transform  works  well  for  separable  factors,  it  does  not 
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always  lead  to  a  stable  digital  realization  starting  with  a 
stable  analog  transfer  function. 

The  following  are  suggested  for  further  research: 

1.  Other  types  of  canonic  factors  need  to  be  considered. 
For  example,  the  general  form: 

H(j»lfjw2)   =   [1  +  (ja31T1)P+(jaJ2T2)q+(jco1T3)r(jw2T4)t]a 

2.  The  work  should  be  extended  to  non-causal  (half- 
plane)  filters. 

3.  Transfer  functions  which  have  circular  symmetry  in 
the  discrete  domain  need  to  be  investigated. 

4.  Another  approach  is  to  consider  half-plane 
recursive  algorithms,  by  working  backwards  from  the  discrete 
domain  (z, ,zj  to  continuous  domain  (s,,s2). 

5.  Stability  of  canonic  factors  needs  to  be  investigated 
in  detail . 

6.  The  use  of  predistortion  should  also  be  considered. 
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